wave_calibration.c

00001 /******************************************************************************
00002 * E.S.O. - VLT project
00003 *
00004 * 
00005 *
00006 * who       when      what
00007 * --------  --------  ----------------------------------------------
00008 * schreib  13/07/00  created
00009 */
00010 
00011 /************************************************************************
00012 *   NAME
00013 *     wave_calibration.c -
00014 *     routines needed for wavelength calibration
00015 *
00016 *   SYNOPSIS
00017 * 
00018 *   1) FitParams ** newFitParams( int n_params )
00019 *
00020 *   2) void destroyFitParams ( FitParams ** params )
00021 *
00022 *   3) void dumpFitParamsToAscii ( FitParams ** params, char * filename )
00023 *
00024 *   4) void dumpAsciiToFitParams ( FitParams ** params, char * filename )
00025 *
00026 *   5) int findLines( OneImage * lineImage, 
00027 *                     float    * wave_position, 
00028 *                     float    * wave_intensity,
00029 *                     int        n_lines, 
00030 *                     int     ** row_clean,
00031 *                     float   ** wavelength_clean,
00032 *                     float      beginWave, 
00033 *                     float      dispersion,
00034 *                     float      mindiff, 
00035 *                     int        halfWidth,
00036 *                     int      * n_found_lines,
00037 *                     float      sigma,
00038 *                     int      * sum_lines )
00039 *
00040 *   6) int readList( char * listname, float * lineCenter, float * lineIntensity )
00041 *
00042 *
00043 *   7) int linefit ( OneImage  *  mergedImage, 
00044 *                    FitParams *  par,
00045 *                    float        fwhm,
00046 *                    int          lineInd,
00047 *                    int          column, 
00048 *                    int          halfWidth, 
00049 *                    int          lineRow,
00050 *                    float        min_amplitude )
00051 *
00052 *   8) int fitLines ( OneImage  *  line_image, 
00053 *                     FitParams ** allParams,
00054 *                     float        fwhm,
00055 *                     int       *  n_lines, 
00056 *                     int       ** row,
00057 *                     float     ** wavelength, 
00058 *                     int          width,
00059 *                     float        min_amplitude ) 
00060 *
00061 *   9) float polyfit( FitParams ** par, 
00062 *                     int          column,
00063 *                     int          n_lines,
00064 *                     int          n_rows,
00065 *                     float        dispersion,
00066 *                     float        max_residual,
00067 *                     float *      acoefs, 
00068 *                     float *      dacoefs, 
00069 *                     int   *      n_reject,
00070 *                     int          n_fitcoefs )
00071 *
00072 *   10) float coefsCrossFit ( int      n_columns,
00073 *                            float *  acoefs,
00074 *                            float *  dacoefs, 
00075 *                            float *  bcoefs,
00076 *                            int      n_fitcoefs,
00077 *                            float    sigma_factor )
00078 *
00079 *
00080 *  11) OneImage * waveMap( OneImage * lineImage,
00081 *                          float   ** bcoefs,
00082 *                          int        n_a_fitcoefs,
00083 *                          int        n_b_fitcoefs,
00084 *                          float    * wavelength,
00085 *                          float    * intensity,
00086 *                          int        n_lines,
00087 *                          int        magFactor,
00088 *                          int      * bad_column_mask,
00089 *                          int        n_bad_columns )
00090 *
00091 *  12) int wavelengthCalibration( OneImage   * image, 
00092 *                                 FitParams ** par ,
00093 *                                 float     ** bcoefs,
00094 *                                 float      * wave,
00095 *                                 int          n_lines,
00096 *                                 int       ** row_clean,
00097 *                                 float     ** wavelength_clean,
00098 *                                 int        * n_found_lines,
00099 *                                 float        dispersion,
00100 *                                 int          halfWidth,
00101 *                                 float        minAmplitude,
00102 *                                 float        max_residual,
00103 *                                 float        fwhm,
00104 *                                 int          n_a_fitcoefs,
00105 *                                 int          n_b_fitcoefs,
00106 *                                 float        sigmaFactor )
00107 *
00108 *  13) OneImage * convolveImageByGauss( OneImage * lineImage,
00109 *                                       int        hw )
00110 *
00111 *  14) OneImage * definedResampling( OneImage * image,
00112 *                                    OneImage * calimage,
00113 *                                    int        n_params,
00114 *                                    int        n_rows,
00115 *                                    double   * dispersion,
00116 *                                    float    * minval,
00117 *                                    float    * maxval,
00118 *                                    double   * centralLambda,
00119 *                                    int      * centralpix )
00120 *
00121 *   DESCRIPTION
00122 *
00123 *   1) allocates memory for a new vector of 
00124 *      FitParams data structures
00125 *   2) frees memory of a vector of FitParams data structures
00126 *   3) dumps the fit parameters to an ASCII file
00127 *   4) dumps ASCII information to an allocated FitParams data structure
00128 *   5) determines the pixel shift between the line list 
00129 *      and the real image by using the beginning wavelength
00130 *      on the detector and the dispersion estimate.
00131 *   6) reads the line data of the calibration lamps
00132 *   7) fits a gaussian to a 1-dimensional slice of an image, 
00133 *      this routine uses the routine lsqfit_c as a non-linear
00134 *      least square fit method (Levenberg-Marquardt).               
00135 *   8) calculates and stores the fit parameters of the neon 
00136 *      emission lines of a neon frame by using the linefit 
00137 *      routine.
00138 *   9) fits a second order polynom 
00139 *      lambda[i] = a1 + a2*pos[i] + a3*pos[i]^2
00140 *      to determine the connection between the listed wave-
00141 *      length values and the gauss-fitted positions for each
00142 *      image column using the singular value decomposition 
00143 *      method. 
00144 *  10) Fits the each single parameter of the three fit parameters 
00145 *      acoefs from polyfit through the image columns
00146 *  11) this routine determines a wavelength calibration map 
00147 *      frame associating a wavelength value to each pixel
00148 *      by using the fit coefficients determined before.
00149 *  12) this routine takes an image from a calibration
00150 *      emission lamp and delivers the fit coefficients of  
00151 *      a polynomial fit across the columns 
00152 *      of the coefficients of the polynomial line position
00153 *      fits as output. Furthermore it delivers an array of the fit parameters
00154 *      as output. This routine expects Nyquist sampled spectra 
00155 *     (either an interleaved image or an image convolved with an 
00156 *      appropriate function in spectral direction)
00157 *  13) convolves an emission line image with a gaussian
00158 *      with user given integer half width by using the eclipse 
00159 *      routine function1d_filter_lowpass().
00160 *  14) Given a source image and a corresponding wavelength
00161 *      calibration file this routine produces an image
00162 *      in which elements in a given row are associated
00163 *      with a single wavelength. It thus corrects for 
00164 *      the wavelength shifts between adjacent elements
00165 *      in the rows of the input image. The output image
00166 *      is larger in the wavelength domain than the input
00167 *      image with pixels in each column corresponding to 
00168 *      undefined (blank, ZERO) values. The distribution
00169 *      of these undefined values varies from column to
00170 *      column. The input image is resampled at discrete
00171 *      wavelength intervals using a polynomial interpolation
00172 *      routine. 
00173 *      The wavelength intervals (dispersion) and the 
00174 *      central wavelength are defined and stable for each
00175 *      used grating. Thus, each row has a defined wavelength
00176 *      for each grating. Only the number of rows can be 
00177 *      changed by the user. 
00178 *
00179 *   FILES
00180 *
00181 *   ENVIRONMENT
00182 *
00183 *   RETURN VALUES 
00184 *
00185 *   CAUTIONS 
00186 *
00187 *   EXAMPLES
00188 *
00189 *   SEE ALSO
00190 *
00191 *   BUGS   
00192 *
00193 *------------------------------------------------------------------------
00194 */
00195 
00196 #include "vltPort.h"
00197 
00198 /* 
00199  * System Headers
00200  */
00201 
00202 /* 
00203  * Local Headers
00204  */
00205 
00206 #include "wave_calibration.h"
00207 #include "solve_poly_root.h"
00208 
00209 
00210 /*----------------------------------------------------------------------------
00211    Function     :       newFitParams()
00212    In           :       number of spectral lines that will be fitted
00213    Out          :       allocated array of FitParams data structures
00214    Job          :       allocates memory for a new array of 
00215                         FitParams data structures
00216  ---------------------------------------------------------------------------*/
00217 
00218 FitParams ** newFitParams( int n_params )
00219 {
00220     FitParams ** new_params =NULL;
00221     FitParams  * temp_params =NULL;
00222     float * temp_fit_mem =NULL;
00223     float * temp_derv_mem=NULL;
00224     int i ;
00225 
00226     
00227     if ( n_params <= 0 )
00228     {
00229         cpl_msg_error ("newFitParams:"," wrong number of lines to fit\n") ;
00230         return NULL ;
00231     }
00232 
00233     if ( NULL == (new_params = (FitParams **) cpl_calloc ( n_params , sizeof (FitParams*) ) ) )
00234     {
00235         cpl_msg_error ("newFitParams:"," could not allocate memory\n") ;
00236         return NULL ;
00237     }
00238     if ( NULL == (temp_params = cpl_calloc ( n_params , sizeof (FitParams) ) ) )
00239     {
00240         cpl_msg_error ("newFitParams:"," could not allocate memory\n") ;
00241         return NULL ;
00242     }
00243     if ( NULL == (temp_fit_mem    = (float *) cpl_calloc( n_params*MAXPAR, sizeof (float) ) ) )
00244     {
00245         cpl_msg_error ("newFitParams:"," could not allocate memory\n") ;
00246         return NULL ;
00247     }
00248 
00249 
00250     if ( NULL == (temp_derv_mem   = (float *) cpl_calloc( n_params*MAXPAR, sizeof (float) ) ) )
00251     {
00252         cpl_msg_error ("newFitParams:"," could not allocate memory\n") ;
00253         return NULL ;
00254     }
00255 
00256     for ( i = 0 ; i < n_params ; i ++ )
00257     {
00258         new_params[i] = temp_params+i;
00259         new_params[i] -> fit_par    = temp_fit_mem+i*MAXPAR;
00260         new_params[i] -> derv_par   = temp_derv_mem+i*MAXPAR;
00261         new_params[i] -> column     = 0 ;
00262         new_params[i] -> line       = 0 ;
00263         new_params[i] -> wavelength = 0. ;
00264         new_params[i] -> n_params   = n_params ;
00265     }
00266 
00267     return new_params ;
00268 }
00269 
00270 /*----------------------------------------------------------------------------
00271    Function     :       destroyFitParams()
00272    In           :       fit params to destroy
00273    Out          :       nothing
00274    Job          :       frees memory of an array of FitParams data structures
00275  ---------------------------------------------------------------------------*/
00276 
00277 void destroyFitParams ( FitParams ** params )
00278 { 
00279     
00280     if ( params == NULL )
00281     {
00282         return ;
00283     }
00284 
00285     cpl_free ( params[0] -> fit_par ) ;
00286     cpl_free ( params[0] -> derv_par ) ;
00287     cpl_free ( params[0] ) ;
00288     cpl_free ( params ) ;
00289 }
00290 
00291 /*----------------------------------------------------------------------------
00292    Function     :       dumpFitParamsToAscii()
00293    In           :       params: fit params to dump
00294                         filename
00295    Out          :       filled ASCII file
00296    Job          :       dumps the fit parameters to an ASCII file
00297  ---------------------------------------------------------------------------*/
00298 
00299 void dumpFitParamsToAscii ( FitParams ** params, char * filename )
00300 {
00301     FILE * fp ;
00302     int    i ;
00303 
00304     if ( NULL == params )
00305     {
00306         cpl_msg_error ("dumpFitParamsToAscii:"," no fit parameters available!\n") ;
00307         return ;
00308     }
00309 
00310     if ( NULL == filename )
00311     {
00312         cpl_msg_error ("dumpFitParamsToAscii:"," no filename available!\n") ;
00313         return ;
00314     }
00315 
00316     if ( NULL == (fp = fopen ( filename, "w" ) ) )
00317     {
00318         cpl_msg_error("dumpFitParamsToAscii:"," cannot open %s\n", filename) ;
00319         return ;
00320     }
00321 
00322     for ( i = 0 ; i < params[0] -> n_params ; i++ )
00323     {
00324         fprintf(fp, "%d %d %d %f %f %f %f %f %f %f %f %f\n", 
00325                      params[i]->n_params, params[i]->column, params[i]->line, 
00326                      params[i]->wavelength, params[i]->fit_par[0], params[i]->fit_par[1], 
00327                      params[i]->fit_par[2], params[i]->fit_par[3],
00328                      params[i]->derv_par[0], params[i]->derv_par[1], 
00329                      params[i]->derv_par[2], params[i]->derv_par[3] ) ; 
00330     }
00331     fclose(fp) ;
00332 }
00333 
00334 /*----------------------------------------------------------------------------
00335    Function     :       dumpAsciiToFitParams()
00336    In           :       allocated dummy for the fit params 
00337                         filename
00338    Out          :       params: filled FitParams object
00339    Job          :       dumps ASCII information to an allocated FitParams data structure
00340  ---------------------------------------------------------------------------*/
00341 
00342 void dumpAsciiToFitParams ( FitParams ** params, char * filename )
00343 {
00344     FILE * fp ;
00345     int i ;
00346 
00347     if ( NULL == params )
00348     {
00349         cpl_msg_error ("dumpAsciiToFitParams:"," no fit parameters available!\n") ;
00350         return ;
00351     }
00352 
00353     if ( NULL == filename )
00354     {
00355         cpl_msg_error ("dumpAsciiToFitParams:"," no filename available!\n") ;
00356         return ;
00357     }
00358 
00359     if ( NULL == (fp = fopen ( filename, "r" ) ) )
00360     {
00361         cpl_msg_error("dumpAsciiToFitParams:"," cannot open %s\n", filename) ;
00362         return ;
00363     }
00364     
00365     for ( i = 0 ; i < params[0]->n_params ; i++ )
00366     {
00367         fscanf(fp, "%d %d %d %f %f %f %f %f %f %f %f %f\n", 
00368                     &params[i]->n_params, &params[i]->column, &params[i]->line, 
00369                     &params[i]->wavelength, &params[i]->fit_par[0], &params[i]->fit_par[1], 
00370                     &params[i]->fit_par[2], &params[i]->fit_par[3],
00371                     &params[i]->derv_par[0], &params[i]->derv_par[1], 
00372                     &params[i]->derv_par[2], &params[i]->derv_par[3] ) ; 
00373     }
00374     fclose(fp) ;
00375 }
00376 
00377 /*----------------------------------------------------------------------------
00378    Function     :       findLines()
00379    In           :       lineImage: merged emission line image,
00380                         wave_position: wavelength list in Angstroems
00381                         wave_intensity: corresponding intensity list 
00382                         n_lines: number of lines in list
00383                         row_clean: resulting list of the row indices but without the 
00384                                    lines that are too close to each other for the fit
00385                         wavelength_clean: corrected wavelength list corresponding to
00386                                           the row_clean array
00387                         beginWave: beginning wavelength on detector in microns 
00388                         dispersion: dispersion of the grating on the detector
00389                                     (microns per pixel, attention: merged image).
00390                         mindiff: minimal difference of mean and median column
00391                                  intensity to do the correlation. 
00392                                  This is done to avoid correlations in columns 
00393                                  without emission line intensity.
00394                         halfWidth:   half width of the box where the line must sit,
00395                         n_found_lines: number of found and correlated 
00396                                        emission lines in a column.
00397                         sigma: sigma of Gaussian that is convolved with the artificial spectrum
00398    Out          :       0 if all went o.k.
00399                         row: resulting list of the row indices of the line positions
00400                         row_clean: resulting list of the row indices but without the 
00401                                    lines that are too close to each other for the fit
00402                         wavelength: wavelength from the list corresponding to the 
00403                                     found row positions
00404                         wavelength_clean: corrected wavelength list corresponding to
00405                                           the row_clean array
00406                         n_found_lines: number of found and correlated 
00407                                        emission lines in a column.
00408                         sum_lines: total sum of found and correlated emission lines
00409                         -1 if something has gone wrong 
00410    Job          :       determines the pixel shift between the line list 
00411                         and the real image by using the beginning wavelength
00412                         on the detector and the dispersion estimate.
00413  ---------------------------------------------------------------------------*/
00414 
00415 int findLines( OneImage * lineImage, 
00416                float    * wave_position, 
00417                float    * wave_intensity,
00418                int        n_lines, 
00419                int     ** row_clean,
00420                float   ** wavelength_clean,
00421                float      beginWave, 
00422                float      dispersion1,
00423                float      dispersion2,
00424                float      mindiff, 
00425                int        halfWidth,
00426                int      * n_found_lines,
00427                float      sigma,
00428                int      * sum_lines )
00429 {
00430     int     ** row ;
00431     float   ** wavelength ;
00432     float buf1, buf2 ;
00433     float meanval ;
00434     float colmedian ;
00435     float * column, * tempcol ;
00436     float * lines ;
00437     float * conv_lines ;
00438     float * wave_buffer ;
00439     float * wave_mem;
00440     int   * dummy_row ;
00441     int i, j, k, m ;
00442     int position ;
00443     int gmax, gmin ;
00444     int col ;
00445     int * row_mem;
00446     float sum ;
00447     float angst ;
00448 
00449     if ( NullImage == lineImage )
00450     {
00451         cpl_msg_error ("findLines:"," no image given\n") ;
00452         return -1 ;
00453     }
00454 
00455     if ( n_lines <= 0 || NULL == wave_position ) 
00456     {
00457         cpl_msg_error("findLines:"," no line list given\n") ;
00458         return -1 ;
00459     }
00460     if ( NULL == wave_intensity ) 
00461     {
00462         cpl_msg_error("findLines:"," no intensity list given\n") ;
00463         return -1 ;
00464     }
00465      
00466     if ( dispersion1 == 0. )
00467     {
00468         cpl_msg_error("findLines:"," wrong dispersion given\n") ;
00469         return -1 ;
00470     }
00471    
00472     if ( row_clean == NULL )
00473     {
00474         cpl_msg_error("findLines:"," row_clean array is not allocated\n") ;
00475         return -1 ;
00476     }
00477 
00478     if ( wavelength_clean == NULL )
00479     {
00480         cpl_msg_error("findLines:"," wavelength_clean array is not allocated\n") ;
00481         return -1 ;
00482     }
00483 
00484     if ( beginWave <= 0. )
00485     {
00486         cpl_msg_error ("findLines:"," impossible beginWave given\n") ;
00487         return -1 ;
00488     }
00489     if ( mindiff < -100. )
00490     {
00491         cpl_msg_error ("findLines:"," wrong mindiff value\n") ;
00492         return -1 ;
00493     }
00494 
00495     if ( halfWidth <= 0 )
00496     {
00497         cpl_msg_error("findLines:"," wrong half width of fit box given\n") ;
00498         return -1 ;
00499     }
00500 
00501     if ( n_found_lines == NULL )
00502     {
00503         cpl_msg_error("findLines:"," n_found_lines not allocated\n") ;
00504         return -1 ;
00505     }
00506 
00507     if ( sigma <= 0. || sigma >= lineImage -> ly / 2)
00508     {
00509         cpl_msg_error("findLines:"," wrong sigma given\n") ;
00510         return -1 ;
00511     }
00512     
00513     /* allocate memory */
00514     row        = (int**) cpl_calloc( lineImage->lx, sizeof(int*)) ;
00515     wavelength = (float**) cpl_calloc( lineImage->lx, sizeof(float*)) ;
00516     row_mem = cpl_calloc( n_lines*lineImage->lx, sizeof(int) ) ;
00517     wave_mem = cpl_calloc( n_lines*lineImage->lx, sizeof(float) ) ;
00518     for ( i = 0 ; i <lineImage -> lx ; i++ )
00519     {
00520         row[i] = row_mem + i*n_lines;
00521         wavelength[i] = wave_mem + i*n_lines;
00522     }
00523 
00524     /* find if the wavelength is given in microns, nanometers or Angstroem */
00525     if ( wave_position[0] > 10000. )
00526     {
00527     /* Angstroem */
00528         angst = 10000. ;
00529     }
00530     else if ( wave_position[0] > 1000. && wave_position[0] < 10000. )
00531     {
00532     /* nanometers */
00533     angst = 1000. ;
00534     }
00535     else
00536     {
00537     /* microns */
00538     angst = 1. ;
00539     }
00540 
00541     /*-------------------------------------------------------------------------------
00542      * compute the mean and median intensity value in the given column
00543      * and determine if there is enough intensity in the column to do the correlation
00544      */
00545     tempcol = (float*) cpl_calloc(lineImage->ly, sizeof(float)) ;
00546     *sum_lines = 0 ;
00547     buf1 = 0. ;
00548     buf2 = 0. ;
00549     /* allocate memory */
00550     column      = (float*) cpl_calloc(lineImage -> ly, sizeof (float)) ;
00551     lines       = (float*) cpl_calloc(lineImage -> ly, sizeof (float)) ;
00552     conv_lines  = (float*) cpl_calloc(lineImage -> ly, sizeof (float)) ;
00553     wave_buffer = (float*) cpl_calloc(lineImage -> ly, sizeof (float)) ;
00554     dummy_row   = (int*)   cpl_calloc(n_lines, sizeof(int)) ;
00555 
00556     for ( col = 0 ; col < lineImage->lx ; col++ )
00557     {
00558         n_found_lines[col] = 0 ;
00559         sum = 0. ;
00560         for ( i = 0 ; i < lineImage -> ly ; i++ )
00561         {
00562             if (qfits_isnan(lineImage->data[col + i*lineImage->lx]) )
00563             {
00564                 tempcol[i] = 0. ;
00565                 continue ;
00566             }
00567 
00568         sum = sum + lineImage -> data[col + i*lineImage -> lx] ;
00569         tempcol[i] = lineImage -> data[col + i*lineImage -> lx];
00570         }
00571         meanval = sum / lineImage -> ly ;
00572         /* lets assume the median is the background */
00573         colmedian =  median ( tempcol, lineImage->ly);
00574 
00575         if ( meanval - colmedian < mindiff )
00576         {
00577         cpl_msg_warning("findLines:"," sorry, not enough intensity (mean: %f, diff: %f) in image column: %d to correlate emission lines\n", meanval, meanval - colmedian,col) ;
00578         continue ;
00579     }
00580 
00581         for ( i = 0 ; i < lineImage -> ly ; i++ )
00582     {
00583       conv_lines[i]=0;
00584       wave_buffer[i]=0;
00585     }
00586     for ( i = 0 ; i < n_lines ; i++ )
00587     {
00588       dummy_row[i] = 0;
00589     }
00590 
00591         /* go through the column with index col */
00592         for ( i = 0 ; i < lineImage -> ly ; i++ )
00593         {
00594             if ( qfits_isnan(lineImage->data[col+i*lineImage->lx]) )
00595             {
00596                 column[i] = 0. ;
00597             }
00598             else
00599             {
00600                 column[i] = lineImage -> data[col + i*lineImage -> lx] ;
00601             }
00602 
00603             /* determine the line position on the pixels */
00604             /*lines[i] = (dispersion * (float) i + beginWave) * angst ;*/ 
00605             lines[i] = (dispersion1 * (float) (i-lineImage->ly/2) + dispersion2 * (float) (i-lineImage->ly/2)  * (float) (i-lineImage->ly/2) + beginWave) * angst ; 
00606 
00607             /* ---------------------------------------------------------------
00608              * find the nearest line position for each wavelength in the list 
00609              * and set this position to the given line intensity as weight 
00610              */
00611             for ( j = 0 ; j < n_lines ; j ++ )
00612             {
00613                 buf1 = 0. ;
00614                 buf2 = 0. ;
00615                 if ( (wave_position[j] >= (lines[i] - fabs(dispersion1)/2.*angst)) && 
00616                      (wave_position[j] <= (lines[i] + fabs(dispersion1)/2.*angst)) ) 
00617                 {
00618                     buf1 = wave_intensity[j] ; /* set the given line intensity as weight */
00619                     buf2 = wave_position[j] ;
00620                     break ;
00621                 }
00622             }
00623             lines[i] = buf1 ;
00624             wave_buffer[i] = buf2 ;      /* get the wavelength associated with the corresponding 
00625                                            found emission line */
00626     
00627             /* convolve the artificial spectrum by a Gaussian with given sigma value */
00628             if ( lines[i] != 0. )
00629             {
00630                 /* consider only +- 2 sigma */
00631                 gmin = nint((float) i - 2. * sigma) ;
00632                 gmax = nint((float) i + 2. * sigma) ;
00633 
00634                 /* be aware of image margins */
00635                 if ( gmin < 0 )
00636                 {
00637                     gmin = 0 ;
00638                 }
00639                 if ( gmax >= lineImage -> ly )
00640                 {
00641                     gmax = lineImage -> ly - 1 ;
00642                 } 
00643                 for ( j = gmin ; j <= gmax ; j++ )
00644                 {
00645                     conv_lines[j] += 
00646                         lines[i] * exp( (double)( -0.5*((j - i)*(j - i)))/(double) (sigma*sigma) ) ; 
00647                 }
00648             }
00649         }
00650 
00651         /* do the cross correlation */
00652         position = INT32_MAX ;
00653         position = correlation(column, conv_lines, lineImage -> ly ) ; 
00654         if ( abs (position) > lineImage -> ly / 4 )
00655         {
00656             cpl_msg_warning("findLines:"," sorry, shift of artificial data relative to image (%d) seems to be too high in column: %d\n", position, col) ;
00657             continue ;
00658         }
00659 
00660         j = 0 ;
00661         for ( i = 0 ; i < lineImage -> ly ; i ++ )
00662         {
00663             if ( lines[i] != 0.0 )
00664             {
00665                 if ( (i - position) >= 0 && (i - position) < lineImage -> ly )
00666                 {
00667                     row[col][j] = i - position ;
00668                     /* get the wavelength corresponding to found line row index */
00669                     wavelength[col][j] = wave_buffer[i] / angst ;
00670                     j++ ;
00671                 }
00672             }
00673         }
00674 
00675 
00676         /* -------------------------------------------------------------------------------------------
00677          *  determine the row_clean array, that means, take only the row values if the distance between
00678          *  adjacent lines is large enough for the fit
00679          */
00680         for ( k = 1 ; k <= j ; k ++ )
00681         {
00682             if (dummy_row[k-1] != -1)
00683             {
00684                 dummy_row[k-1] = row[col][k-1] ;
00685             }
00686             if ( (row[col][k] - row[col][k-1]) <= 2*halfWidth )
00687             {
00688                 dummy_row[k-1] = -1 ;
00689         if (k<n_lines)
00690                   dummy_row[k]   = -1 ;
00691             }
00692             if ( (row[col][k+1] - row[col][k]) <= 2*halfWidth)
00693             {
00694         if (k<n_lines)
00695                   dummy_row[k]   = -1 ;
00696         if (k+1<n_lines)
00697                   dummy_row[k+1] = -1 ;
00698             }
00699         }
00700 
00701         m = 0 ;
00702         for ( k = 0 ; k < j ; k ++ )
00703         {
00704             if ( dummy_row[k] != -1 && dummy_row[k] != 0 )
00705             {
00706                 row_clean[col][m] = dummy_row[k] ;
00707                 wavelength_clean[col][m] = wavelength[col][k] ;
00708                 m ++ ;
00709             }
00710         }
00711 
00712         n_found_lines[col] = m ;
00713 
00714         *sum_lines += n_found_lines[col] ;
00715     }
00716     cpl_free (column) ;
00717     cpl_free (lines) ;
00718     cpl_free (conv_lines) ;
00719     cpl_free (dummy_row) ;
00720     cpl_free (wave_buffer) ;
00721     cpl_free (row_mem) ;
00722     cpl_free (wave_mem) ;
00723     cpl_free (tempcol) ;
00724     cpl_free (row) ;
00725     cpl_free (wavelength) ;
00726 
00727     return 0 ;
00728 }
00729 
00730 /*----------------------------------------------------------------------------
00731    Function     :       readList()
00732    In           :       name of the list file, arrays to store the wavelength
00733                         and the intensities of the emission lines
00734    Out          :       number of lines in file
00735    Job          :       reads the line data of the calibration lamps
00736  ---------------------------------------------------------------------------*/
00737 
00738 int readList( char * listname, float * lineCenter, float * lineIntensity )
00739 {
00740     FILE * fp ;
00741     int i, n_lines ;
00742     
00743     if ( NULL == lineCenter )
00744     {
00745         cpl_msg_error("readList:"," lineCenter array is not allocated\n") ;
00746         return -1 ;
00747     }
00748 
00749     if ( NULL == lineIntensity )
00750     {
00751         cpl_msg_error("readList:"," lineIntensity array is not allocated\n") ;
00752         return -1 ;
00753     }
00754 
00755     if ( NULL == (fp = fopen ( listname, "r" ) ) )
00756     {
00757         cpl_msg_error("readList:"," cannot open %s\n", listname) ;
00758         return -1 ;
00759     }
00760 
00761     i = 0 ;
00762     while ( fscanf( fp, "%f %f", &lineCenter[i], &lineIntensity[i] ) != EOF )
00763     {
00764         i ++ ;
00765     }
00766     n_lines = i ;
00767     fclose(fp) ;
00768            
00769     return n_lines ;
00770 }
00771 
00772 
00773 /*----------------------------------------------------------------------------
00774    Function     :       linefit()
00775    In           :       mergedImage: image of a calibration emission lamp,
00776                         par:         dummys for the resulting fitting parameters,
00777                         fwhm:        guess value for full width of half maximum of gaussian
00778                         lineInd:     index of the emission line,
00779                         column:      present index of the image column,
00780                         halfWidth:   half width of the box where the line must sit,
00781                         lineRow:     row index where the line is expected,
00782                         min_amplitude: minimum line amplitude with respect to the background 
00783                        to do the fit
00784    Out          :       the fitting parameter data structure containing the 
00785                         resulting parameters.
00786                         integers: 
00787                         number of iterations if all was ok,
00788                         -8   if no input image was given,
00789                         -9   if no dummy for the fit parameters is given,
00790                         -10  if the wrong column index was given,
00791                         -11  if the wrong box width was given,
00792                         -12  if the wrong row index was given,
00793                         -13  if a wrong minimum amplitude factor was given
00794                         -14  if the spectral vector data structure memory
00795                              could not be allocated          
00796                         -15  wrong row index or box width was given,
00797                         -16  signal too low to fit,         
00798                         -17  least squares fit failed
00799    Job          :       fits a gaussian to a 1-dimensional slice of an image, 
00800                         this routine uses the routine lsqfit_c as a non-linear
00801                         least square fit method (Levenberg-Marquardt).               
00802  ---------------------------------------------------------------------------*/
00803 
00804 int linefit ( OneImage  *  mergedImage, 
00805               FitParams *  par,
00806               float        fwhm,
00807               int          lineInd,
00808               int          column, 
00809               int          halfWidth, 
00810               int          lineRow,
00811               float        min_amplitude,
00812           Vector    *  line,
00813           int       *  mpar,
00814           float     *  xdat,
00815           float     *  wdat )
00816 {
00817     int i, j ;
00818     int iters, xdim, ndat ;
00819     int numpar, its ;
00820     int position ;
00821     float maxval, tol, lab ;
00822 
00823     if ( mergedImage == NullImage )
00824     {
00825         cpl_msg_error ("linefit:"," no image given as input\n") ;
00826         return -8 ;
00827     }
00828     if ( par == NULL )
00829     {
00830         cpl_msg_error("linefit:"," fit parameters not given\n") ;
00831         return -9 ;
00832     }
00833     if ( column < 0 || column > mergedImage -> lx )
00834     {
00835         cpl_msg_error ("linefit:"," wrong column number\n") ;
00836         return -10 ;
00837     }
00838     if ( halfWidth < 0 || halfWidth > mergedImage -> ly )
00839     {
00840         cpl_msg_error ("linefit:"," wrong width given\n") ;
00841         return -11 ;
00842     }
00843     if ( lineRow < 0 || lineRow > mergedImage -> ly )
00844     {
00845         cpl_msg_error ("linefit:"," wrong number of row of the line given\n") ;
00846         return -12 ;
00847     }
00848     if ( min_amplitude < 1. )
00849     {
00850         cpl_msg_error ("linefit:"," wrong minimum amplitude\n") ;
00851         return -13 ;
00852     }
00853     
00854     /* initialise the Vector */
00855     for ( i = 0 ; i < line -> n_elements ; i++) 
00856     {
00857         line->data[i] = 0;
00858     }
00859     
00860     par -> column = column  ;
00861     par -> line   = lineInd ;
00862 
00863     /* determine the values of the spectral vector given as input */
00864     /* go through the chosen column */
00865 
00866     j = 0 ;
00867     for ( i = lineRow-halfWidth ; i <= lineRow+halfWidth ; i++ ) 
00868     {
00869         if ( i < 0 || i >= mergedImage -> ly )
00870         {
00871             cpl_msg_error ("linefit:"," wrong line position or width given\n") ;
00872             return -15 ;
00873         }
00874         else
00875         {
00876             line -> data[j] = mergedImage -> data[column + i*mergedImage -> lx] ;
00877             j ++ ;
00878         }
00879     } 
00880 
00881     /*-------------------------------------------------------------------- 
00882      * go through the spectral vector 
00883      * determine the maximum pixel value in the spectral vector 
00884      */
00885     maxval = -FLT_MAX ;
00886     position = -INT32_MAX ;
00887     for ( i = 0 ; i < line -> n_elements ; i++ )
00888     {
00889         xdat[i] = i ;
00890         wdat[i] = 1.0 ;
00891         if ( line -> data[i] >= maxval )
00892         {
00893             maxval = line -> data[i] ;
00894             position = i ;
00895         }
00896     }
00897 
00898     /* set initial values for the fitting routine */
00899     xdim     = XDIM ;
00900     ndat     = line -> n_elements ;
00901     numpar   = MAXPAR ;
00902     tol      = TOL ;
00903     lab      = LAB ;
00904     its      = ITS ;
00905     par -> fit_par[1] = fwhm ;
00906     par -> fit_par[2] = (float) position ;
00907     par -> fit_par[3] = (float) (line -> data[0] + line -> data[line->n_elements - 1]) / 2.0 ;
00908     par -> fit_par[0]  = maxval - (par -> fit_par[3]) ;
00909 
00910     /* exclude low signal cases */
00911     if ( par->fit_par[0] < min_amplitude )
00912     {
00913         cpl_msg_debug ("linefit:"," sorry, amplitude of line too low to fit: %f\n",par->fit_par[0] ) ;
00914         return -16 ;
00915     }
00916 
00917     for ( i = 0 ; i < MAXPAR ; i++ )
00918     {
00919         par -> derv_par[i] = 0.0 ;
00920         mpar[i] = 1 ;
00921     }
00922 
00923     /* finally, do the least square fit using a gaussian */
00924     if ( 0 > ( iters = lsqfit_c( xdat, &xdim, line -> data, wdat, &ndat, par -> fit_par,  
00925                                  par -> derv_par, mpar, &numpar, &tol, &its, &lab )) )  
00926     {
00927         cpl_msg_debug ("linefit:"," lsqfit_c: least squares fit failed, error no.: %d\n", iters) ; 
00928         return -17 ;
00929     }
00930 
00931     /* correct the fitted position for the given row of the line in image coordinates */
00932     par -> fit_par[2] =  (float) (lineRow - halfWidth) + par -> fit_par[2] ;
00933 
00934     /* all was o.k. */
00935     return iters ;
00936 }
00937 
00938 /*----------------------------------------------------------------------------
00939    Function     :       fitLines
00940    In           :       line_image: merged image of a calibration lamp ,
00941                         allParams:  allocated vector of FitParams data structures,
00942                         fwhm:       guess value for full width of half maximum of gaussian
00943                         n_lines:    number of neon lines that will be fitted in one column ,
00944                         row:        list of the rows of the fitted lines
00945                         wavelength: list of wavelength corresponding to the found line rows
00946                         width:      half width of a box around the found rows within the line is fitted
00947                         min_amplitude: minimum line amplitude with respect to the background 
00948                        to do the fit
00949    Out          :       filled FitParams data structure vector,
00950                         number of successfully fitted lines,
00951                         errors: negative integers resulting from the linefit
00952                                 routine and:
00953                         -18: no image given,
00954                         -19: number of emission lines or number of slitlets is wrong,
00955                         -20: vector of the slitlet boundaries or of the line rows or
00956                              of the half width are empty.
00957                         -21: no wavelength array given.
00958    Job          :       calculates and stores the fit parameters of the neon 
00959                         emission lines of a neon frame by using the linefit 
00960                         routine.
00961  ---------------------------------------------------------------------------*/
00962 
00963 int fitLines ( OneImage  *  line_image, 
00964                FitParams ** allParams,
00965                float        fwhm,
00966                int       *  n_lines, 
00967                int       ** row,
00968                float     ** wavelength, 
00969                int          width,
00970                float        min_amplitude ) 
00971 {
00972     int i, k, l ;
00973     int result ;
00974     Vector * line;
00975     int    * mpar;
00976     float  * xdat, * wdat;
00977 
00978     if ( line_image == NULL )
00979     {
00980         cpl_msg_error ("fitLines:"," no image given\n") ;
00981         return -18 ;
00982     }
00983     if ( n_lines == NULL )
00984     {
00985         cpl_msg_error ("fitLines:"," no counter of emission lines\n") ;
00986         return -19 ;
00987     } 
00988     if ( row == NULL || width <= 0 )
00989     {
00990         cpl_msg_error ("fitLines:"," row or width vectors are empty\n") ;
00991         return -20 ;
00992     }
00993     if ( wavelength == NULL )
00994     {
00995         cpl_msg_error ("fitLines:"," no wavelength array given\n") ;
00996         return -21 ;
00997     }
00998 
00999     k = 0 ;
01000 
01001     /* allocate memory for the spectral vector */
01002     line = newVector (2*width + 1) ;
01003     /* allocate memory */
01004     xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
01005     wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
01006     mpar = (int *)   cpl_calloc( MAXPAR, sizeof (int) ) ;
01007     
01008     /* go through the columns */
01009     for ( i = 0 ; i < line_image -> lx ; i++ )
01010     {
01011         if ( n_lines[i] == 0 )
01012         {
01013             continue ;
01014         }
01015         /* go through the emission lines in a column */
01016         for ( l = 0 ; l < n_lines[i] ; l++ )
01017         {
01018             if ( row[i][l] <= 0 )
01019             {
01020                 continue ;
01021             }
01022 
01023             /* --------------------------------------------------------------
01024              * fit the single lines using linefit and store the parameters in
01025              * an array of the FitParams data structure allParams[].
01026              */ 
01027             if ( (result = linefit ( line_image, allParams[k], fwhm, l, i, 
01028                                      width, row[i][l], min_amplitude,line,mpar,xdat,wdat ) ) < 0 )
01029             {
01030                 cpl_msg_debug ("fitLines:"," linefit failed, error no.: %d, column: %d, row: %d, line: %d\n", result, i, row[i][l], l) ;
01031                 continue ;
01032             }
01033             if ( (allParams[k] -> fit_par[0] <= 0.) || (allParams[k] -> fit_par[1] <= 0.)
01034                   || (allParams[k] -> fit_par[2] <= 0.) )
01035             {
01036                 cpl_msg_warning ("fitLines:"," negative fit parameters in column: %d, line: %d\n", i, l) ;
01037                 continue ;
01038             }
01039             allParams[k] -> wavelength = wavelength[i][l] ;
01040             k++ ;
01041         }
01042     }
01043 
01044     /* free memory */
01045     destroyVector(line);
01046     cpl_free(xdat);
01047     cpl_free(wdat);
01048     cpl_free(mpar);
01049     
01050     /* all is o.k. */
01051     return k ;
01052 }    
01053 
01054 /*----------------------------------------------------------------------------
01055    Function     :       polyfit()
01056    In           :       par:          filled array of fit parameter structure
01057                         column:       image column index
01058                         n_lines:      number of found lines in column
01059                         n_rows:       number of image rows
01060                         dispersion:   microns per pixel
01061                         max_residual: maximum residual value, beyond that value
01062                                       the fit is rejected.
01063                         acoefs:       array of the 3 coefficients of the fitted
01064                                       parabola
01065                         dacoefs:      array of standard deviations of the 3 coefficients
01066                         n_reject:     rejected number of fits due to high residuals
01067                         n_fitcoefs:   number of polynomial coefficients to fit
01068    Out          :       chisq, the three fit coefficients acoefs[i] and their
01069                         standard deviations dacoefs[i], the rejected number of fits
01070                         due to too high residuals: n_reject
01071                         FLT_MAX in case of error
01072    Job          :       fits a polynomial
01073                         lambda[i] = a1 + a2*pos[i] + a3*pos[i]^2 +...
01074                         to determine the connection between the listed wave-
01075                         length values and the gauss-fitted positions for each
01076                         image column using the singular value decomposition
01077                         method.
01078  ---------------------------------------------------------------------------*/
01079 
01080 float polyfit( FitParams ** par,
01081                int          column,
01082                int          n_lines,
01083                int          n_rows,
01084                float        dispersion,
01085                float        max_residual,
01086                float *      acoefs,
01087                float *      dacoefs,
01088                int   *      n_reject,
01089                int          n_fitcoefs )
01090 {
01091     float ** ucoefs, ** vcoefs, ** covar ;
01092     float *mem;
01093     float wcoefs[n_fitcoefs] ;
01094     float * lambda, * posit ;
01095     float * weight, * resid ;
01096     float * newlam, * newpos, * newwet ;
01097     float chisq, result ;
01098     float offset ;
01099     int num, found ;
01100     int i, j, k, n ;
01101 
01102     /* reset the fit coefficients and their errors */
01103     for ( i = 0 ; i < n_fitcoefs ; i++ )
01104     {
01105         acoefs[i]  = 0. ;
01106         dacoefs[i] = 0. ;
01107     }
01108     if ( NULL == par )
01109     {
01110         cpl_msg_error("polyfit:"," no fit params given\n");
01111         return FLT_MAX ;
01112     }
01113 
01114     if ( 0 >= n_lines )
01115     {
01116       /*
01117         cpl_msg_warning ("polyfit:"," sorry, number of lines is wrong") ;
01118       */
01119         return FLT_MAX ;
01120     }
01121     if ( 0 >= n_rows )
01122     {
01123         cpl_msg_error ("polyfit:"," sorry, number of rows is wrong") ;
01124         return FLT_MAX ;
01125     }
01126     if ( dispersion == 0. )
01127     {
01128         cpl_msg_error ("polyfit:"," sorry, wrong dispersion given") ;
01129         return FLT_MAX ;
01130     }
01131 
01132     offset = (float)(n_rows - 1)/2. ;
01133 
01134     /* allocate memory */
01135     
01136     mem = (float*) cpl_calloc( n_lines*7, sizeof (float) ) ;
01137     lambda = mem;
01138     posit  = mem + n_lines;
01139     weight = mem + n_lines*2;
01140     resid  = mem + n_lines*3;
01141     newlam = mem + n_lines*4;
01142     newpos = mem + n_lines*5;
01143     newwet = mem + n_lines*6;
01144     
01145     /*lambda = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01146     posit  = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01147     weight = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01148     resid  = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01149     newlam = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01150     newpos = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01151     newwet = (float*) cpl_calloc( n_lines, sizeof (float) ) ;*/
01152 
01153     /* allocate coefficient matrices*/
01154     ucoefs = matrix ( 1, n_lines, 1, n_fitcoefs ) ;
01155     vcoefs = matrix ( 1, n_lines, 1, n_fitcoefs ) ;
01156     covar  = matrix ( 1, n_fitcoefs, 1, n_fitcoefs ) ;
01157 
01158     /* go through all fit parameters */
01159     n = 0 ;
01160     for ( i = 0 ; i < (par[0] -> n_params) ; i ++ )
01161     {
01162         found = -1 ;
01163         /* find the given column and go through the lines in that column */
01164         for ( j = 0 ; j < n_lines ; j ++ )
01165         {
01166             if ( (par[i] -> column == column) && (par[i] -> line == j) )
01167             {
01168                 found = i ;
01169             }
01170             else
01171             {
01172                 continue ;
01173             }
01174 
01175             /* store only fit params with reasonable values */
01176             if ( par[found] -> derv_par[2] != 0. && par[found] -> fit_par[2] > 0. &&
01177                  par[found] -> wavelength > 0. && par[found] -> fit_par[1] > 0. &&
01178                  par[found] -> fit_par[0] > 0. )
01179             {
01180                 /* ---------------------------------------------------------------------------
01181                  * store the found position, error of the position as weight and the associated
01182                  * wavelength values of the fitted lines
01183                  */
01184                 posit[n]  = par[found] -> fit_par[2] ;
01185                 weight[n] = par[found] -> derv_par[2] ;
01186                 lambda[n] = par[found] -> wavelength ;
01187                 n ++ ;
01188             }
01189             else
01190             {
01191                 continue ;
01192             }
01193         }
01194 
01195     }
01196 
01197     num = n ;
01198     if ( num < n_fitcoefs )
01199     {
01200         cpl_msg_warning("polyfit:"," not enough lines found in column %d to determine the three coefficients.\n", column) ;
01201         for ( i = 0 ; i < n_fitcoefs ; i++ )
01202         {
01203             acoefs[i]  = ZERO ;
01204             dacoefs[i] = ZERO ;
01205         }
01206         free_matrix ( ucoefs, 1/*, n_lines*/,    1/*, n_fitcoefs*/ ) ;
01207         free_matrix ( vcoefs, 1/*, n_lines*/,    1/*, n_fitcoefs*/ ) ;
01208         free_matrix ( covar,  1/*, n_fitcoefs*/, 1/*, n_fitcoefs*/ ) ;
01209         /*cpl_free (lambda) ;
01210         cpl_free (posit) ;
01211         cpl_free (weight) ;
01212         cpl_free (resid) ;
01213         cpl_free (newlam) ;
01214         cpl_free (newpos) ;
01215         cpl_free (newwet) ;*/
01216     cpl_free (mem);
01217         return FLT_MAX ;
01218     }
01219 
01220     /*--------------------------------------------------------------------------------
01221      * scale the pixel position values to smaller than 1 and transform the weights to
01222      * wavelength units 
01223      */
01224 
01225     for ( i = 0 ; i < num ; i ++ )
01226     {
01227         posit[i] = (posit[i] - offset)/offset ;
01228         weight[i] *= fabs(dispersion) ;
01229     }
01230 
01231     /* do the fit using the singular value decomposition method */
01232     svd_fitting( posit - 1, lambda - 1, weight - 1, num, acoefs-1, n_fitcoefs,
01233             ucoefs, vcoefs, wcoefs-1, covar, &chisq, fpol ) ;
01234 
01235     /* scale the linear and the quadratic coefficient */
01236     for ( i = 1 ; i < n_fitcoefs ; i++ )
01237     {
01238         acoefs[i] /= pow(offset, i) ;
01239     }
01240 
01241     /* now that we have determined the fit coefficients, find the residuals */
01242     *n_reject = 0 ;
01243 
01244     j = 0 ;
01245     for ( i = 0 ; i < num ; i++ )
01246     {
01247         result = 0. ;
01248         for ( k = 0 ; k < n_fitcoefs ; k++ )
01249         {
01250             result += acoefs[k] * pow(posit[i], k) ;
01251         }
01252 
01253         resid[i] = lambda[i] - result ;
01254 
01255         if ( fabs( resid[i] ) > max_residual)
01256         {
01257             (*n_reject) ++ ;
01258         }
01259         else
01260         {
01261             newlam[j] = lambda[i] ;
01262             newpos[j] = posit[i] ;
01263             newwet[j] = weight[i] ;
01264             j++ ;
01265         }
01266     }
01267 
01268     num = j ;
01269     if ( num >= n_fitcoefs )
01270     {
01271         svd_fitting( newpos - 1, newlam - 1, newwet - 1, num, acoefs-1, n_fitcoefs, ucoefs,
01272                 vcoefs, wcoefs-1, covar, &chisq, fpol ) ;
01273 
01274         /* scale the resulting coefficients */
01275         for ( i = 0 ; i < n_fitcoefs ; i++ )
01276         {
01277             acoefs[i] /= pow(offset, i) ;
01278             dacoefs[i] = sqrt( (double) covar[i+1][i+1] ) / pow(offset, i) ;
01279         }
01280     }
01281     else
01282     {
01283         cpl_msg_warning ("polyfit:"," too many lines rejected (number: %d) due to high residuals, fit coefficients are set zero, in column: %d\n", *n_reject, column) ;
01284         for ( i = 0 ; i < n_fitcoefs ; i++ )
01285         {
01286             acoefs[i]  = ZERO ;
01287             dacoefs[i] = ZERO ;
01288         }
01289     }
01290 
01291     free_matrix ( ucoefs, 1/*, n_lines*/,    1/*, n_fitcoefs*/ ) ;
01292     free_matrix ( vcoefs, 1/*, n_lines*/,    1/*, n_fitcoefs*/ ) ;
01293     free_matrix ( covar,  1/*, n_fitcoefs*/, 1/*, n_fitcoefs*/ ) ;
01294     /*cpl_free (lambda) ;
01295     cpl_free (posit) ;
01296     cpl_free (weight) ;
01297     cpl_free (resid) ;
01298     cpl_free (newlam) ;
01299     cpl_free (newpos) ;
01300     cpl_free (newwet) ;*/
01301     cpl_free (mem);
01302 
01303     return chisq ;
01304 }
01305 
01306 /*----------------------------------------------------------------------------
01307    Function     :       coefsCrossFit()
01308    In           :       n_columns:    number of image columns
01309                         acoefs:       coeffs fitted in polyfit
01310                                       note: this is a vector of coefficients with
01311                                             the same index for all columns
01312                         dacoefs:      fit errors of the corresponding acoefs
01313                         bcoefs:       the fitted coefs
01314                         n_fitcoefs:   number of fit coefficients
01315                         sigma_factor: factor of sigma beyond which the
01316                                       column coefficients are discarded for the fit
01317    Out          :       chisq, the found fit coefficients
01318    Job          :       Fits each single polnomial coefficient acoefs resulting from polyfit
01319                         across the image columns
01320  ---------------------------------------------------------------------------*/
01321 
01322 float coefsCrossFit ( int      n_columns,
01323                       float *  acoefs,
01324                       float *  dacoefs,
01325                       float *  bcoefs,
01326                       int      n_fitcoefs,
01327                       float    sigma_factor )
01328 {
01329     float col_index, sub_col_index[n_columns] ;
01330     float sub_acoefs[n_columns], sub_dacoefs[n_columns] ;
01331     float wcoefs[n_fitcoefs] ;
01332     float ** ucoefs, **vcoefs, **covar ;
01333     float chisq ;
01334     float * acoefsclean ;
01335     double sum, sumq, mean ;
01336     double sigma ;
01337     double cliphi, cliplo ;
01338     float offset ;
01339     int i, n, num, ndata ;
01340     int nc ;
01341 
01342     if ( n_columns < 1 )
01343     {
01344         cpl_msg_error("coefsCrossFit:"," wrong number of image columns given\n") ;
01345         return FLT_MAX ;
01346     }
01347     if ( acoefs == NULL || dacoefs == NULL )
01348     {
01349         cpl_msg_error("coefsCrossFit:"," coeffs or errors of coefficients are not given\n") ;
01350         return FLT_MAX ;
01351     }
01352     if ( bcoefs == NULL )
01353     {
01354         cpl_msg_error("coefsCrossFit:"," coeffs are not allocated\n") ;
01355         return FLT_MAX ;
01356     }
01357 
01358     if ( n_fitcoefs < 1 )
01359     {
01360         cpl_msg_error("coefsCrossFit:"," wrong number of fit coefficients\n") ;
01361         return FLT_MAX ;
01362     }
01363     if ( sigma_factor <= 0. )
01364     {
01365         cpl_msg_error("coefsCrossFit:"," impossible sigma_factor given!\n") ;
01366         return FLT_MAX ;
01367     }
01368 
01369     offset = (float)(n_columns - 1) / 2. ;
01370 
01371     /* ----------------------------------------------------------
01372      * determine the clean mean and sigma value of the coefficients,
01373      * that means reject 10 % of the extreme low and high values
01374      */
01375     nc = 0 ;
01376     for ( i = 0 ; i < n_columns ; i++ )
01377     {
01378         if ( qfits_isnan(acoefs[i]) || acoefs[i] == 0. || dacoefs[i] == 0. )
01379         {
01380             continue ;
01381         }
01382         else
01383         {
01384             nc++ ;
01385         }
01386     }
01387     acoefsclean = (float*) cpl_calloc(nc , sizeof(float)) ;
01388     nc = 0 ;
01389     for ( i = 0 ; i < n_columns ; i++ )
01390     {
01391         if ( qfits_isnan(acoefs[i]) || acoefs[i] == 0. || dacoefs[i] == 0. )
01392         {
01393             continue ;
01394         }
01395         else
01396         {
01397             acoefsclean[nc] = acoefs[i] ;
01398             nc++ ;
01399         }
01400     }
01401     pixel_qsort(acoefsclean, nc) ;
01402     sum   = 0. ;
01403     sumq  = 0. ;
01404     mean  = 0. ;
01405     sigma = 0. ;
01406     n     = 0 ;
01407     for ( i = (int)((float)nc*LOW_REJECT) ; i < (int)((float)nc*HIGH_REJECT) ; i++ )
01408     {
01409         sum  += (double)acoefsclean[i] ;
01410         sumq += ((double)acoefsclean[i] * (double)acoefsclean[i]) ;
01411         n ++ ;
01412     }
01413     mean          = sum/(double)n ;
01414     sigma         = sqrt( sumq/(double)n - (mean * mean) ) ;
01415     cliphi        = mean + sigma * (double)sigma_factor ;
01416     cliplo        = mean - sigma * (double)sigma_factor ;
01417 
01418     /* fit only the reasonnable values */
01419     num = 0 ;
01420     for ( i = 0 ; i < n_columns ; i++ )
01421     {
01422         /* associate the column indices to the corresponding array */
01423         col_index = (float) i ;
01424 
01425         /* take only the reasonnable coefficients */
01426         if ( !qfits_isnan(acoefs[i]) && (acoefs[i] <= cliphi) && (acoefs[i] >= cliplo) &&
01427              (dacoefs[i] != 0. ) && (acoefs[i] != 0.) )
01428         {
01429             sub_acoefs[num]    = acoefs[i] ;
01430             sub_dacoefs[num]   = dacoefs[i] ;
01431             sub_col_index[num] = col_index ;
01432             num ++ ;
01433         }
01434     }
01435     ndata = num ;
01436 
01437     if ( ndata < n_fitcoefs )
01438     {
01439         cpl_msg_error("coefsCrossFit:"," not enough data found to determine the fit coefficients.\n") ;
01440 
01441         return FLT_MAX ;
01442     }
01443 
01444     /* allocate coefficient matrices */
01445     ucoefs = matrix(1, ndata, 1, n_fitcoefs) ;
01446     vcoefs = matrix(1, ndata, 1, n_fitcoefs) ;
01447     covar  = matrix ( 1, n_fitcoefs, 1, n_fitcoefs ) ;
01448 
01449     /* scale the x-values for the fit */
01450     for ( i = 0 ; i < ndata ; i++ )
01451     {
01452         sub_col_index[i] = (sub_col_index[i] - offset) / offset ;
01453     }
01454 
01455     /* finally, do the singular value decomposition fit */
01456     svd_fitting ( sub_col_index-1, sub_acoefs-1, sub_dacoefs-1, ndata, bcoefs-1,
01457              n_fitcoefs, ucoefs, vcoefs, wcoefs-1, covar, &chisq, fpol ) ;
01458 
01459     /* scale the found coefficients */
01460     for ( i = 0 ; i < n_fitcoefs ; i ++ )
01461     {
01462         bcoefs[i] /= pow(offset, i) ;
01463     }
01464 
01465     /* free memory */
01466     cpl_free (acoefsclean) ;
01467     free_matrix( ucoefs, 1/*, ndata*/, 1/*, n_fitcoefs */) ;
01468     free_matrix( vcoefs, 1/*, ndata*/, 1/*, n_fitcoefs */) ;
01469     free_matrix ( covar, 1/*, n_fitcoefs*/, 1/*, n_fitcoefs*/ ) ;
01470 
01471     return chisq ;
01472 }
01473 
01474 
01475 /*----------------------------------------------------------------------------
01476    Function     :       waveMap()
01477    In           :       lineImage:    image from a calibration emission lamp,
01478                         bcoefs:       transformed fit coefficients
01479                         n_a_fitcoefs: number of fit coefficients for the single 
01480                                       column fits lambda-position
01481                         n_b_fitcoefs: number of fit coefficients for the fits of
01482                                       the single a coefficients across the columns
01483                         wavelength:   wavelength list from lamp file
01484                         intensity:    corresponding line intensity from line list
01485                         n_lines:      number of lines in the list
01486                         magFactor:    magnifying factor of the image for FFT
01487    Out          :       wavelength calibration map image.
01488    Job          :       this routine determines a wavelength calibration map 
01489                         frame associating a wavelength value to each pixel
01490                         by using the fit coefficients determined before.
01491  ---------------------------------------------------------------------------*/
01492 
01493 OneImage * waveMap( OneImage * lineImage,
01494                     float   ** bcoefs,
01495                     int        n_a_fitcoefs,
01496                     int        n_b_fitcoefs,
01497                     float    * wavelength,
01498                     float    * intensity,
01499                     int        n_lines,
01500                     int        magFactor)
01501 {
01502     OneImage * retImage ;
01503     float cenpos, cenpix ;
01504     float centreval, centrepix, wavelag ;
01505     float  pixvalue ;
01506     double a[n_a_fitcoefs],  wave[n_lines] ;
01507     float a_initial ;
01508     int i, j, k, l/*, m*/, line, col, row, found, sign ;
01509     int var, maxlag, cmin, cmax, offset ;
01510     float emline[2*magFactor*lineImage -> ly], spec[2*magFactor*lineImage -> ly] ;
01511     double * result ;
01512     float col_off ;
01513     float angst ;
01514     double xcorr_max ;
01515     int delta ;
01516     double z[2*(n_a_fitcoefs - 1)] ;
01517     gsl_poly_complex_workspace * w ;
01518     
01519     if ( NullImage == lineImage )
01520     {
01521         cpl_msg_error("waveMap:","no image given\n") ;
01522         return NullImage ;
01523     }
01524 
01525     if ( NULL == wavelength || n_lines <= 0 )
01526     {
01527         cpl_msg_error("waveMap:","no wavelength list given\n") ;
01528         return NullImage ;
01529     }
01530 
01531     if ( NULL == intensity )
01532     {
01533         cpl_msg_error("waveMap:","no intensity list given\n") ;
01534         return NullImage ;
01535     }
01536 
01537     if ( NULL == bcoefs )
01538     {
01539         cpl_msg_error("waveMap:","no coefficients given\n") ;
01540         return NullImage ;
01541     }
01542 
01543     if ( magFactor <= 1 )
01544     {
01545         cpl_msg_error("waveMap:","wrong magnifying factor given\n") ;
01546         return NullImage ;
01547     }
01548     
01549     /* allocate memory */
01550     if ( NULL == ( retImage = new_image( lineImage -> lx, lineImage -> ly ) ))
01551     {
01552         cpl_msg_error("waveMap:","cannot allocate a new image\n");
01553         return NullImage ;
01554     }
01555 
01556     var    = (magFactor - 1)*(magFactor - 1) ;
01557     offset = lineImage -> ly * (magFactor/4 + 1) ;
01558 
01559     /* find out if Angstroem or microns are used */
01560     if ( wavelength[0] > 10000. )
01561     {
01562     /* Angstroem */
01563         angst = 10000. ;
01564     }
01565     else if ( wavelength[0] > 1000. && wavelength[0] < 10000. )
01566     {
01567     /* nanometers */
01568     angst = 1000. ;
01569     }
01570     else
01571     {
01572     /* microns */
01573     angst = 1. ;
01574     }
01575 
01576     /* go through the image columns */
01577     for ( col = 0 ; col < lineImage -> lx ; col++ )
01578     {
01579         /* initialize the emline array for each column */
01580         for ( i = 0 ; i < 2*magFactor*lineImage -> ly ; i++ )
01581         {
01582             emline[i] = 0. ;
01583         }
01584         col_off = (float)col - (float)(lineImage->lx-1)/2. ;
01585 
01586         /* determine the coefficients by using the given bcoefs */
01587         for ( i = 0 ; i < n_a_fitcoefs ; i++ ) 
01588         {
01589             /* initialize coefficients and solution */
01590             a[i] = 0. ;
01591             if (i < n_a_fitcoefs-1)
01592             {
01593                 z[2*i] = 0. ;
01594                 z[2*i+1] = 0. ;
01595             }
01596             for ( j = 0 ; j < n_b_fitcoefs ; j++ )
01597             {
01598                 a[i] += bcoefs[i][j] * pow(col_off, j) ;
01599             }
01600         }
01601         a_initial = a[0] ;
01602         
01603         /* go through the lines and generate an artificial spectrum */
01604         for ( line = 0 ; line < n_lines ; line++ )
01605         {
01606             /* go from Angstroem to micron */
01607             wave[line] = wavelength[line]/angst ;
01608 
01609             /* ---------------------------------------------------------------------
01610              * solve the polynomial for the exact offset of the line that means
01611              * find the root of the polynomial of order n_fitcoefs - 1
01612              */
01613             a[0] = a_initial - wave[line] ;
01614 
01615             if (NULL == (w = gsl_poly_complex_workspace_alloc(n_a_fitcoefs)))
01616             {
01617                 cpl_msg_error("waveMap:","could not allocate complex workspace!") ;
01618                 destroy_image(retImage) ;
01619                 return NullImage ;
01620             }
01621             if (-1 == gsl_poly_complex_solve(a, n_a_fitcoefs, w, z))
01622             {
01623                 cpl_msg_error("waveMap:","gsl_poly_complex_solve did not work!") ;
01624                 destroy_image(retImage) ;
01625                 return NullImage ;
01626             }
01627             gsl_poly_complex_workspace_free(w) ;
01628 
01629             
01630             j = 0 ;
01631             found = -1 ;
01632             for ( i = 0 ; i < n_a_fitcoefs - 1  ; i++ )
01633             {
01634                 /* test for appropriate solution */
01635                 if( z[2*i] > (-1.)*(float) lineImage->ly/2. && 
01636                     z[2*i] < (float)lineImage -> ly/2. && z[2*i+1] == 0. )
01637                 {
01638                     found = 2*i ;
01639                     j ++ ;
01640                 }
01641                 else
01642                 {
01643                     continue ;
01644                 } 
01645             }
01646             if ( j == 0 )
01647             {
01648                 cpl_msg_warning("waveMap:","no offset solution found for line %d in column %d\n", line, col) ;
01649                 continue ;
01650             } 
01651             else if ( j == 1 )
01652             {
01653                 cenpos = z[found] + (float) lineImage->ly /2. ;
01654             }
01655             else
01656             {
01657                 cpl_msg_warning("waveMap:","two or more offset solutions found for line %d in column %d\n", line, col) ;
01658                 continue ;
01659             }
01660              
01661             /*---------------------------------------------------------------------------------
01662              * magnify image by the given factor add an additional offset 
01663              */
01664             cenpix = cenpos * (float) magFactor + (float) offset ;  
01665             
01666             /* determine max and min pixel limits over which line should be convolved */
01667             cmin = (nint(cenpix) - (var-1)) > 0 ? nint(cenpix) - (var-1) : 0 ;
01668             cmax = (nint(cenpix) + (var-1)) < 2*magFactor*lineImage -> ly ? 
01669                     nint(cenpix) + (var-1) : 2*magFactor*lineImage -> ly ;
01670 
01671             /* convolve neon lines with Gaussian function */
01672             for ( j = cmin ; j < cmax ; j++ )
01673             {
01674                 emline[j] += intensity[line] * exp((double)(-0.5*(j-cenpix)*(j-cenpix))/(double)var) ;
01675             }
01676         }
01677                 
01678         /*----------------------------------------------------------------------------------------- 
01679          * for each column, map the image data points onto an magFactor times bigger element grid for FFT
01680          * in the cross correlation, first initialize the two helping arrays for each new column.
01681          */
01682         for ( k = 0 ; k < 2*magFactor*lineImage -> ly ; k++ )
01683         {
01684             spec[k] =  0. ;
01685         }
01686    
01687         /* now take the image data points of the column and put them into the spec array */
01688         for ( row = 0 ; row < lineImage -> ly ; row++ ) /* go through the column */
01689         {
01690             /* insert 8 values for each image row (magnification) and add same offset as for emline array */
01691             for ( l = 0 ; l < magFactor ; l++ )   
01692         {
01693             /* set bad pixels or negative values to zero */
01694                 if (!qfits_isnan(lineImage -> data[col + row*lineImage -> lx]) &&
01695                     (lineImage -> data[col + row*lineImage -> lx] > 0.))
01696                 {
01697                     spec[offset + l + (row * magFactor)] = 
01698                 lineImage -> data[col + row*lineImage -> lx] ;     
01699                 }
01700                 else
01701                 {
01702                     spec[offset + l + (row * magFactor)] = 0. ;
01703                 }
01704             }
01705     }
01706  
01707         /* now call the cross correlation routine */
01708         if (NULL == (result = xcorrel(spec, 2*magFactor*lineImage->ly, emline, 2*magFactor*lineImage -> ly, 
01709                          magFactor*lineImage->ly, &delta, &maxlag, &xcorr_max)) ) 
01710         {
01711         cpl_msg_warning ("cross-correl","cross correlation did not work, col: %d is set ZERO\n", col) ;
01712             for ( row = 0 ; row < lineImage -> ly ; row++ )
01713             {
01714                 retImage -> data[col + row*lineImage -> lx] = ZERO ;
01715             }
01716             continue ;
01717         }
01718     
01719         if ( xcorr_max <= 0. )
01720         {
01721             cpl_msg_warning ("cross-correl","cross correlation sum is negative, col: %d is set ZERO\n", col) ;
01722             for ( row = 0 ; row < lineImage -> ly ; row++ )
01723             {
01724                 retImage -> data[col + row*lineImage -> lx] = ZERO ;
01725             }
01726             cpl_free(result) ;
01727             continue ;
01728         }
01729 
01730         wavelag = (float) -delta / (float) magFactor ;
01731         if ( fabs(wavelag) > (float)lineImage->ly/20. )
01732         {
01733             cpl_msg_warning ("cross-correl","wave lag too big, col: %d is set ZERO\n", col) ;
01734             for ( row = 0 ; row < lineImage -> ly ; row++ )
01735             {
01736                 retImage -> data[col + row*lineImage -> lx] = ZERO ;
01737             }
01738             cpl_free(result) ;
01739         continue ;
01740         }
01741 
01742         /*----------------------------------------------------------------------------- 
01743          * determine new zero order coefficient centreval, of which the formula is
01744          * determined by setting equal a polynomial shifted by wavelag with the same 
01745          * higher order coefficients and set the new zero order coefficient to 
01746          * get both sides of the equation approximately equal.
01747          */ 
01748         centreval = a_initial ;
01749         for ( i = 1 ; i < n_a_fitcoefs ; i++ )
01750         {
01751             if ( i%2 == 0 )
01752             {
01753                 sign = -1 ;
01754             }
01755             else
01756             {
01757                 sign = 1 ;
01758             }
01759             centreval += (float)sign * a[i]*pow(wavelag, i) ;
01760         }
01761 
01762         /* prepare to write out wavelength as pixel values */
01763         for ( row = 0 ; row < lineImage -> ly ; row++ )
01764         {
01765             centrepix = (float)row - ((float)lineImage->ly - 1.)/2. ;
01766             pixvalue = 0. ;
01767             for ( i = 1 ; i < n_a_fitcoefs ; i++ )
01768             {
01769                 pixvalue += a[i]*pow(centrepix, i) ; 
01770             }
01771             retImage -> data[col + row*lineImage -> lx] = 
01772                 centreval + pixvalue ; 
01773         }
01774         cpl_free(result) ;   
01775     }
01776     return retImage ;
01777 }                    
01778 
01779 /*----------------------------------------------------------------------------
01780    Function     :       wavelengthCalibration()
01781    In           :       image:        merged image from a calibration emission lamp,
01782                         wave:         wavelength array read from the wavelength list
01783                         n_lines:      number of lines in the wavelength list
01784                         row_clean: resulting list of the row indices but without the
01785                                    lines that are too close to each other for the fit
01786                                    output of findLines()
01787                         wavelength_clean: corrected wavelength list corresponding to
01788                                           the row_clean array
01789                                           output of findLines()
01790                         n_found_lines: output of findLines(): total number of found emission lines
01791                         dispersion:   dispersion of spectrum: micron per pixel
01792                         halfWidth:    half width of the box where the line must sit
01793                         minAmplitude: minimum amplitude of the Gaussian to do the fit
01794                         max_residual: maximum residual value, beyond that value
01795                                       the polynomial lambda-position fit is rejected.
01796                         fwhm:         first guess for the full width of half maximum
01797                                       of the gaussian line fit
01798                         n_a_fitcoefs: number of fit coefficients for the single
01799                                       column fits: lambda-position
01800                         n_b_fitcoefs: number of fit coefficients for the fits of
01801                                       the single a coefficients across the columns
01802                         sigmaFactor:  factor of the standard deviation of the determined
01803                                       polynomial coefficients of the columns beyond
01804                                       which these coefficients are not used to carry out
01805                                       the polynomial fit across the columns.
01806                         pixel_tolerance: maximum tolerated difference between estimated
01807                                          and fitted line positions.
01808    Out          :       0 if all went o.k., -1 if something went wrong.
01809                         bcoefs: array of cooefficients of the polynomial fit across the columns.
01810                         par: array of the resulting FitParams data structure
01811    Job          :       this routine takes an image from a calibration
01812                         emission lamp and delivers the fit coefficients of
01813                         a polynomial fit across the columns
01814                         of the coefficients of the polynomial line position
01815                         fits as output. Furthermore it delivers an array of the fit parameters
01816                         as output. This routine expects Nyquist sampled spectra
01817                         (either an interleaved image or an image convolved with an
01818                          appropriate function in spectral direction)
01819  ---------------------------------------------------------------------------*/
01820 
01821 int wavelengthCalibration( OneImage   * image,
01822                            FitParams ** par ,
01823                            float     ** bcoefs,
01824                            float      * wave,
01825                            int          n_lines,
01826                            int       ** row_clean,
01827                            float     ** wavelength_clean,
01828                            int        * n_found_lines,
01829                            float        dispersion,
01830                            int          halfWidth,
01831                            float        minAmplitude,
01832                            float        max_residual,
01833                            float        fwhm,
01834                            int          n_a_fitcoefs,
01835                            int          n_b_fitcoefs,
01836                            float        sigmaFactor,
01837                            float        pixel_tolerance )
01838 
01839 {
01840     int          i, j, k ;
01841     int          n_fit ;
01842     int          n_reject ;
01843     float     *  acoefs ;
01844     float     *  dacoefs ;
01845     float     ** abuf ;
01846     float     ** dabuf ;
01847     float        chisq_poly, chisq_cross ;
01848     int          zeroind ;
01849     /*float     *  mem ;*/
01850 
01851 
01852     if (  NullImage == image )
01853     {
01854         cpl_msg_error("wavelengthCalibration:","no image given\n") ;
01855         return -1 ;
01856     }
01857     if ( par == NULL )
01858     {
01859         cpl_msg_error("wavelengthCalibration:","no fit parameter data structure given\n") ;
01860         return -1 ;
01861     }
01862     if ( wave == NULL )
01863     {
01864         cpl_msg_error("wavelengthCalibration:","no wavelength list given\n") ;
01865         return -1 ;
01866     }
01867     if ( n_lines <= 0 )
01868     {
01869         cpl_msg_error("wavelengthCalibration:","impossible number of lines in line list given\n") ;
01870         return -1 ;
01871     }
01872     if ( row_clean == NULL )
01873     {
01874         cpl_msg_error("wavelengthCalibration:","no row_clean array given\n") ;
01875         return -1 ;
01876     }
01877     if ( wavelength_clean == NULL )
01878     {
01879         cpl_msg_error("wavelengthCalibration:","no wavelength_clean array given\n") ;
01880         return -1 ;
01881     }
01882 
01883     if ( dispersion == 0. )
01884     {
01885         cpl_msg_error("wavelengthCalibration:","impossible dispersion given\n") ;
01886         return -1 ;
01887     }
01888 
01889     if ( halfWidth <= 0 || halfWidth > image->ly/2 )
01890     {
01891         cpl_msg_error("wavelengthCalibration:","impossible half width of the fitting box given\n") ;
01892         return -1 ;
01893     }
01894     if ( minAmplitude < 1. )
01895     {
01896         cpl_msg_error("wavelengthCalibration:","impossible minimal amplitude\n") ;
01897         return -1 ;
01898     }
01899 
01900     if ( max_residual <= 0. || max_residual > 1. )
01901     {
01902         cpl_msg_error("wavelengthCalibration:","impossible max_residual given\n") ;
01903         return -1 ;
01904     }
01905 
01906     if ( fwhm <= 0. || fwhm > 10. )
01907     {
01908         cpl_msg_error("wavelengthCalibration:","impossible fwhm given\n") ;
01909 
01910         return -1 ;
01911     }
01912 
01913     if ( n_a_fitcoefs <= 0 || n_a_fitcoefs > 9 )
01914     {
01915         cpl_msg_error("wavelengthCalibration:","unrealistic n_a_fitcoefs given\n") ;
01916         return -1 ;
01917     }
01918 
01919     if ( n_b_fitcoefs <= 0 || n_b_fitcoefs > 9 )
01920     {
01921         cpl_msg_error("wavelengthCalibration:"," unrealistic n_b_fitcoefs given\n") ;
01922         return -1 ;
01923     }
01924     if ( sigmaFactor <= 0. )
01925     {
01926         cpl_msg_error("wavelengthCalibration:"," impossible sigmaFactor given\n") ;
01927         return -1 ;
01928     }
01929 
01930     /* initialize the variables */
01931     n_reject = 0 ;
01932     n_fit = 0 ;
01933 
01934     /* fit each found line by using a gaussian function and determine the exact position */
01935     if ( 0 > (n_fit = fitLines( image , par, fwhm, n_found_lines, row_clean, wavelength_clean,
01936                                 halfWidth, minAmplitude )) )
01937     {
01938         cpl_msg_error("wavelengthCalibration:"," cannot fit the lines, error code of fitLines: %d\n", n_fit) ;
01939         return -1 ;
01940     }
01941 
01942     /* first check for faked lines like bad pixels */
01943     if ( -1 == checkForFakeLines (par, dispersion, wavelength_clean, row_clean, n_found_lines,
01944                                   image->lx, pixel_tolerance) )
01945     {
01946         cpl_msg_error("waveCal:"," cannot fit the lines, error code of fitLines: %d\n", n_fit) ;
01947         return -1 ;
01948     }
01949 
01950     /* allocate memory */
01951     if ( NULL == (acoefs = (float*) cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
01952          NULL == (dacoefs = (float*) cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
01953          NULL == (abuf = (float**) cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
01954          NULL == (dabuf = (float**) cpl_calloc (n_a_fitcoefs, sizeof(float*))) )
01955     {
01956         cpl_msg_error("wavelengthCalibration:"," cannot allocate memory\n") ;
01957         return -1 ;
01958     }
01959 
01960     for ( i = 0 ; i < n_a_fitcoefs ; i++ )
01961     {
01962         if ( NULL == (abuf[i] = (float*) cpl_calloc(image -> lx, sizeof(float))) ||
01963              NULL == (dabuf[i] = (float*) cpl_calloc(image -> lx, sizeof(float))) )
01964         {
01965             cpl_msg_error("wavelengthCalibration:"," cannot allocate memory\n") ;
01966             cpl_free(abuf) ;
01967             cpl_free(dabuf) ;
01968             return -1 ;
01969         }
01970     }
01971 
01972     /* fit wavelengths to the corresponding found positions for each column */
01973     k = 0 ;
01974     
01975     for ( i = 0 ; i < image -> lx ; i++ )
01976     {
01977         zeroind = 0 ;
01978         if ( FLT_MAX == (chisq_poly = polyfit( par, i, n_found_lines[i], image -> ly, dispersion,
01979                                                max_residual, acoefs, dacoefs, &n_reject, n_a_fitcoefs)) )
01980         {
01981       /* 
01982            cpl_msg_warning ("wavelengthCalibration:"," error in polyfitt in column: %d\n", i) ;
01983        */
01984             for ( j = 0 ; j < n_a_fitcoefs ; j++ )
01985             {
01986                 acoefs[j] = ZERO ;
01987                 dacoefs[j] = ZERO ;
01988             }
01989         }
01990 
01991         for ( j = 0 ; j < n_a_fitcoefs ; j++ )
01992         {
01993             if ( acoefs[0] <= 0. || acoefs[1] ==0. ||
01994                  dacoefs[j] == 0. || qfits_isnan(acoefs[j]) )
01995             {
01996                 zeroind = 1 ;
01997 
01998             }
01999         }
02000         for ( j = 0 ; j < n_a_fitcoefs ; j++ )
02001         {
02002             if ( zeroind == 0 )
02003             {
02004                 abuf[j][i]  = acoefs[j] ;
02005                 dabuf[j][i] = dacoefs[j] ;
02006             }
02007             else
02008             {
02009                 abuf[j][i]  = ZERO ;
02010                 dabuf[j][i] = ZERO ;
02011             }
02012         }
02013     }
02014 
02015     /* fit each acoefs across the columns to smooth the result */
02016     for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02017     {
02018         if ( FLT_MAX == (chisq_cross = coefsCrossFit( image -> lx, abuf[i], dabuf[i],
02019                                                       bcoefs[i], n_b_fitcoefs, sigmaFactor )) )
02020         {
02021             cpl_msg_error ("wavelengthCalibration:"," cannot carry out the fitting of coefficients across the columns, for the coefficient with index: %d\n", i) ;
02022             for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02023             {
02024                 cpl_free (abuf[i]) ;
02025                 cpl_free (dabuf[i]) ;
02026             }
02027             cpl_free ( acoefs ) ;
02028             cpl_free ( dacoefs ) ;
02029             cpl_free ( abuf ) ;
02030             cpl_free ( dabuf ) ;
02031             return -1 ;
02032         }
02033     }
02034 
02035     /* free all allocated memory */
02036     for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02037     {
02038         cpl_free (abuf[i]) ;
02039         cpl_free (dabuf[i]) ;
02040     }
02041     cpl_free ( acoefs ) ;
02042     cpl_free ( dacoefs ) ;
02043     cpl_free ( abuf ) ;
02044     cpl_free ( dabuf ) ;
02045 
02046     return 0 ;   
02047 }
02048 
02049                            
02050 /*----------------------------------------------------------------------------
02051    Function     :       convolveImageByGauss()
02052    In           :       lineImage:  emission line image  
02053                         hw:         kernel half width of the gaussian 
02054                                     response function
02055    Out          :       emission line image convolved with a gaussian
02056    Job          :       convolves an emission line image with a gaussian
02057                         with user given integer half width by using the eclipse 
02058                         routine function1d_filter_lowpass().
02059  ---------------------------------------------------------------------------*/
02060 
02061 OneImage * convolveImageByGauss( OneImage * lineImage,
02062                                  int        hw )
02063 {
02064     OneImage * returnImage ;
02065     float column_buffer[lineImage->ly] ;
02066     float * filter ;
02067     int col, row ;
02068     
02069     if ( lineImage == NullImage )
02070     {
02071         cpl_msg_error("convolveImageByGauss:"," no input image given!\n") ;
02072         return NullImage ;
02073     }
02074 
02075     if ( hw < 1 )
02076     {
02077         cpl_msg_error("convolveImageByGauss:"," wrong half width given!\n") ;
02078         return NullImage ;
02079     }
02080 
02081     /* allocate memory for returned image */
02082     if ( NULL == ( returnImage = new_image( lineImage -> lx, lineImage -> ly ) ))
02083     {
02084         cpl_msg_error("convolveImageByGauss:"," cannot allocate a new image\n");
02085         return NullImage ;
02086     }
02087 
02088     /* go through the image columns and save them in a buffer */
02089     for ( col = 0 ; col < lineImage -> lx ; col++ )
02090     { 
02091         for ( row = 0 ; row < lineImage -> ly ; row++ )
02092         {
02093             column_buffer[row] = lineImage -> data[col + row*lineImage->lx] ;
02094         }
02095          
02096         /*--------------------------------------------------------------------- 
02097          * now low pass filter the columns by the gaussian and fill the return
02098          * image.
02099          */  
02100         filter = function1d_filter_lowpass( column_buffer,
02101                                             lineImage -> ly,
02102                                             LOW_PASS_GAUSSIAN,
02103                                             hw ) ;
02104         for ( row = 0 ; row < lineImage -> ly ; row++ )
02105         {
02106             returnImage -> data[col + row*lineImage->lx] = filter[row] ;
02107         }
02108         function1d_del(filter) ;
02109     }
02110                         
02111     return returnImage ;
02112 }
02113 
02114 /*----------------------------------------------------------------------------
02115    Function     :       definedResampling()
02116    In           :       image:      source image to be calibrated
02117                         calimage:   wavelength map image 
02118                         n_params:   number of fit parameters for the polynomial interpolation 
02119                                     standard should be 3
02120                                     that means order of polynom + 1
02121                         n_rows:     desired number of rows for the final image,
02122                                     this will be the final number of spectral pixels
02123                                     in the final data cube.
02124                         dispersion: dummy for the resulting dispersion
02125                         minval:     dummy for minimal wavelength value,
02126                         maxval:     dummy for maximal wavelength value 
02127                         centralLambda: dummy for the final central wavelength
02128    Out          :       wavelength calibrated source image, 
02129                         dispersion: resulting spectral dispersion (microns/pixel)
02130                                     is chosen as the minimum dispersion found in 
02131                                     the wavelength map - 2% of this value
02132                         minval:     minimal wavelength value,
02133                         maxval:     maximal wavelength value 
02134                         centralLambda: final central wavelength value
02135                         centralpix:    row of central wavelength 
02136                                        (in image coordinates!)
02137    Job          :       Given a source image and a corresponding wavelength
02138                         calibration file this routine produces an image
02139                         in which elements in a given row are associated
02140                         with a single wavelength. It thus corrects for 
02141                         the wavelength shifts between adjacent elements
02142                         in the rows of the input image. The output image
02143                         is larger in the wavelength domain than the input
02144                         image with pixels in each column corresponding to 
02145                         undefined (blank, ZERO) values. The distribution
02146                         of these undefined values varies from column to
02147                         column. The input image is resampled at discrete
02148                         wavelength intervals using a polynomial interpolation
02149                         routine. 
02150                         The wavelength intervals (dispersion) and the 
02151                         central wavelength are defined and stable for each
02152                         used grating. Thus, each row has a defined wavelength
02153                         for each grating. Only the number of rows can be 
02154                         changed by the user. 
02155  ---------------------------------------------------------------------------*/
02156 
02157 OneImage * definedResampling( OneImage * image,
02158                               OneImage * calimage,
02159                               int        n_params,
02160                               int*       n_rows,
02161                               double   * dispersion,
02162                               float    * minval,
02163                               float    * maxval,
02164                               double   * centralLambda,
02165                               int    * centralpix )
02166 {
02167     OneImage * retImage ;
02168     OneImage * tempCalImage ;
02169     OneImage * tempImage ;
02170     float lambda ;
02171     float dif, lambda_renorm ;
02172     float imagecol[image->ly] ;
02173     float retimagecol[2560] ; /* retimagecol[n_rows] ; */
02174     float calcol[calimage->ly] ;
02175     float x_renorm[n_params] ;
02176     float * imageptr ;
02177     float sum, new_sum ;
02178     float disp, mindisp ;
02179     int calcolpos[2560];
02180     int i/*, j*/, col, row, testrow ;
02181     int half_width, firstpos ;
02182     int dispInd ;
02183     int n ;
02184     int flag;
02185     float temprow;
02186     float minLambda = 0. ;
02187     /*dpoint list[n_params] ;*/
02188     /*double * polycoeffs ;*/
02189     double poly ;
02190     /*float error;*/
02191     int zeroind ;
02192 
02193 
02194     if ( NULL == image )
02195     {
02196         cpl_msg_error("definedResampling:"," source image not given\n") ;
02197         return NullImage ;
02198     }
02199 
02200     if ( NULL == calimage )
02201     {
02202         cpl_msg_error("definedResampling:"," wavelength map image not given\n") ;
02203         return NullImage ;
02204     }
02205 
02206     if ( image -> lx != calimage -> lx ||
02207          image -> ly != calimage -> ly )
02208     {
02209         cpl_msg_error("definedResampling:"," source image and wavelength map image are not compatible in size\n") ;
02210         return NullImage ;
02211     }                              
02212   
02213     if ( n_params < 1 )
02214     {
02215         cpl_msg_error ("definedResampling:"," wrong number of fit parameters given\n") ;
02216         return NullImage ;
02217     }
02218 
02219     if ( n_params > 4 )
02220     {
02221         cpl_msg_warning("definedResampling:"," attention: very high number of fit parameters given, not tested !!!\n") ;
02222     }
02223 
02224     /*if ( n_rows <= calimage -> ly)
02225     {
02226         cpl_msg_error ("definedResampling:"," number of rows of resampled image will be smaller than in wavelength calibration map, information would get lost!\n") ;
02227         return NullImage ;
02228     }*/
02229    
02230     dispInd = 0 ;
02231     /* first determine the dispersion direction */
02232     for ( col = 0 ; col < calimage->lx ; col++ )
02233     {
02234         if ( qfits_isnan(calimage -> data[col]) || calimage -> data[col] <= 0. )
02235         {
02236             continue ;
02237         }
02238         if ((calimage->data[col] - calimage->data[col+(calimage->lx)*(calimage->ly-1)]) > 0. )
02239         {
02240            dispInd--  ;
02241         }
02242         else if ((calimage->data[col] - calimage->data[col+(calimage->lx)*(calimage->ly-1)]) < 0. )
02243         {
02244            dispInd++ ;
02245         }
02246         else
02247         {
02248             continue ;
02249         }
02250     }
02251     if ( dispInd == 0 )
02252     {
02253         cpl_msg_error("definedResampling:"," zero dispersion?\n");
02254         return NullImage ;
02255     }
02256    
02257     /* mirror the wavelength map and the raw image if the dispersion is negative */
02258     if ( dispInd < 0 )
02259     {
02260         /* allocate a temp image */
02261         if ( NULL == ( tempCalImage = new_image( calimage->lx, calimage->ly ) ))
02262         {
02263             cpl_msg_error("definedResampling:"," cannot allocate a new image\n");
02264             return NullImage ;
02265         }
02266         if ( NULL == ( tempImage = new_image( image->lx, image->ly ) ))
02267         {
02268             cpl_msg_error("definedResampling:"," cannot allocate a new image\n");
02269             destroy_image(tempCalImage) ;
02270             return NullImage ;
02271         }
02272         for ( col = 0 ; col < calimage->lx ; col++ )
02273         {
02274             n = calimage->ly - 1 ;
02275             for ( row = 0 ; row < calimage->ly ; row++ )
02276             {
02277                 tempCalImage->data[col+row*calimage->lx] = calimage->data[col+n*calimage->lx] ;
02278                 tempImage->data[col+row*calimage->lx] = image->data[col+n*calimage->lx] ;
02279                 n-- ;
02280             }
02281         }
02282         for ( i = 0 ; i < (int) image->nbpix ; i++ )
02283         {
02284             image->data[i] = tempImage->data[i] ;
02285             calimage->data[i] = tempCalImage->data[i] ;
02286         }
02287         destroy_image(tempCalImage) ;
02288         destroy_image(tempImage) ;
02289     }
02290    
02291     /* determine the max and min pixel value in the first and the last row */
02292     *maxval = -FLT_MAX ;
02293     *minval =  FLT_MAX ;
02294     mindisp = FLT_MAX ;
02295     for ( col = 0 ; col < calimage->lx ; col++ )
02296     {
02297         if ( qfits_isnan(calimage -> data[col]) || calimage -> data[col] <= 0. )
02298         {
02299            continue ;
02300         }
02301         disp = (calimage->data[col+(calimage->lx)*((calimage->ly)-1)]
02302                - calimage->data[col]) / (float)calimage->ly ;
02303         if ( mindisp > disp )
02304         {
02305             mindisp = disp ;
02306         }
02307         if ( *minval >= calimage -> data[col] )
02308         {
02309             *minval = calimage -> data[col] ;
02310         }
02311         if ( *maxval <= calimage -> data[col + (calimage->lx)*((calimage->ly)-1)] )
02312         {
02313             *maxval = calimage -> data[col + (calimage->lx)*((calimage->ly)-1)] ;
02314         }
02315     }
02316     /* find the used grating and set the dispersion to the defined value */
02317     if (*minval > 1.9 )
02318     {
02319         if ( calimage -> ly > 1024 && calimage -> ly < 3000)
02320         {
02321             *dispersion = DISPERSION_K_DITH ;
02322             *centralLambda = CENTRALLAMBDA_K ;
02323         }
02324         else if ( calimage -> ly < 2000)
02325         {
02326             *dispersion = DISPERSION_K ;
02327             *centralLambda = CENTRALLAMBDA_K ;
02328         }
02329         else
02330         {
02331             *dispersion = DISPERSION_K_DITH/2 ;
02332             *centralLambda = CENTRALLAMBDA_K ;
02333         }
02334     }
02335     else if (*minval < 1.2 )
02336     {
02337         if ( calimage -> ly > 1024 )
02338         {
02339             *dispersion = DISPERSION_J_DITH ;
02340             *centralLambda = CENTRALLAMBDA_J ;
02341         }
02342         else
02343         {
02344             *dispersion = DISPERSION_J ;
02345             *centralLambda = CENTRALLAMBDA_J ;
02346         }
02347     }
02348     else 
02349     {
02350         if ( *maxval > 2.3 )
02351         {
02352             if ( calimage -> ly > 1024 )
02353             {
02354                 *dispersion = DISPERSION_HK_DITH ;
02355                 *centralLambda = CENTRALLAMBDA_HK ;
02356             }
02357             else
02358             {
02359                 *dispersion = DISPERSION_HK ;
02360                 *centralLambda = CENTRALLAMBDA_HK ;
02361             }
02362         }
02363         else 
02364         {
02365             if ( calimage -> ly > 1024 )
02366             {
02367                 *dispersion = DISPERSION_H_DITH ;
02368                 *centralLambda = CENTRALLAMBDA_H ;
02369             }
02370             else
02371             {
02372                 *dispersion = DISPERSION_H ;
02373                 *centralLambda = CENTRALLAMBDA_H ;
02374             }
02375         }
02376     }
02377     /*if ( *minval + (float)n_rows * *dispersion < *maxval ) 
02378     {
02379         cpl_msg_error("definedResampling:"," given number of rows too small!\n");
02380         return NullImage ;
02381     }*/
02382     if ( (*maxval - *minval) / *dispersion < (float)calimage->ly ) 
02383     {
02384         cpl_msg_error("definedResampling:"," must be something wrong with the wavelength map!\n");
02385         return NullImage ;
02386     }
02387    
02388     /* determine the central pixel and the lambda in the first image row */
02389     *n_rows = floor(floor(0.5+(*maxval - *minval) / *dispersion)/2+0.5)*2;
02390     *centralpix = *n_rows / 2 ;
02391     minLambda  = *centralLambda - *dispersion * (float)*centralpix ;
02392     /*if ( (minLambda + *dispersion * n_rows) < *maxval ) 
02393     {
02394         cpl_msg_error("definedResampling:"," not enough rows defined \n");
02395         return NullImage ;
02396     }
02397     if ( minLambda  > *minval ) 
02398     {
02399         cpl_msg_error("definedResampling:"," not enough rows defined \n");
02400         return NullImage ;
02401     }*/
02402 
02403     /* allocate memory */
02404     if ( NULL == ( retImage = new_image( image -> lx, *n_rows ) ))
02405     {
02406         cpl_msg_error("definedResampling:"," cannot allocate a new image\n");
02407         return NullImage ;
02408     }
02409     
02410     /* now go through the columns */
02411     for ( col = 0 ; col < retImage -> lx ; col++ )
02412     {
02413         
02414         /*------------------------------------------------------------------ 
02415          * copy the columns of the source image and the wavemap image into
02416          * buffer arrays to speed things up
02417          */
02418         sum = 0. ;
02419         for ( row = 0 ; row < image -> ly ; row++ )
02420         {
02421             imagecol[row] = image    -> data[col + row*image->lx] ; 
02422             if (!qfits_isnan(imagecol[row]))
02423             {
02424                 sum += imagecol[row] ;
02425             }
02426             calcol[row]   = calimage -> data[col + row*calimage->lx] ; 
02427         }
02428         for ( row = 0 ; row < retImage -> ly ; row++ )
02429         {
02430             retimagecol[row] = 0. ;
02431         calcolpos[row] = -1;
02432         }
02433     
02434     for ( row=0 ; row < calimage->ly ; row++)
02435     {
02436         temprow = (calcol[row]- minLambda)/ *dispersion;
02437         if (temprow >= 0 && temprow < retImage->ly)
02438             calcolpos[(int) temprow]  = row;
02439     }
02440         
02441     zeroind = 0 ;
02442        
02443 
02444 
02445         for ( row = 0 ; row < retImage -> ly ; row++ )
02446         {
02447             lambda = minLambda + *dispersion * (float) row ;
02448       
02449             /*--------------------------------------------------------------- 
02450              * lambda must lie between the two available wavelength extremes
02451              * otherwise the image pixels are set to ZERO 
02452              */
02453             if ( row < calimage->ly )
02454             {
02455                 if ( qfits_isnan(calcol[row]) )
02456                 {
02457                     zeroind = 1 ;
02458                 } 
02459             }
02460 
02461             if ( (lambda < calcol[0]) || (lambda > calcol[(calimage->ly)-1]) || zeroind == 1 )
02462             {
02463                 retimagecol[row] = ZERO ;
02464                 continue ;
02465             }
02466 
02467             /*testrow = 0 ; 
02468             while ( lambda > calcol[testrow] )
02469             {
02470                 testrow++ ;
02471             }*/
02472         if (calcolpos[row]==-1) {
02473                 if(row>= (*n_rows-1)) calcolpos[row] = calcolpos[row-1];
02474                 if(row<  (*n_rows-1)) calcolpos[row] = calcolpos[row+1];
02475         }
02476         if (lambda-calcol[calcolpos[row]-1]==0.) calcolpos[row]=calcolpos[row]-1;
02477         testrow = calcolpos[row];
02478 
02479             /*--------------------------------------------------------------------
02480              * at this point calcol[testrow-1] < lambda <= calcol[testrow] 
02481              * now determine the box position in which the polint fit is carried through.
02482              * the half width of the box is half the number of fit parameters.
02483              * Now we determine the start position of the fitting box and treat
02484              * the special case of being near the edge.
02485              */
02486 
02487             if ( n_params % 2 == 0 )
02488             {
02489                 half_width = (int)(n_params/2) - 1 ;
02490             }
02491             else
02492             {
02493                 half_width = (int)(n_params/2) ;
02494             }
02495 
02496 
02497             if ( qfits_isnan(imagecol[testrow]) )
02498             {
02499                 for ( i = row-half_width ; i < row-half_width+n_params ; i++ )
02500                 { 
02501                     if (i < 0) continue ;
02502                     if ( i >= retImage->ly ) continue  ;
02503                     retimagecol[i] = ZERO ;
02504                 }
02505                 imagecol[testrow] = 0. ;
02506             }
02507 
02508         }
02509        
02510 
02511 
02512         /* now loop over the rows and establish the lambda for each row */
02513         new_sum = 0. ;
02514         for ( row = 0 ; row < retImage -> ly ; row++ )
02515         {
02516             if ( qfits_isnan(retimagecol[row]) )
02517             {
02518                 continue ;
02519             }
02520         lambda = minLambda + *dispersion * (float) row ;
02521       
02522             /*--------------------------------------------------------------- 
02523              * lambda must lie between the two available wavelength extremes
02524              * otherwise the image pixels are set to ZERO 
02525              */
02526             if ( (lambda < calcol[0]) || (lambda > calcol[(calimage->ly)-1]) ) 
02527             {
02528                 retimagecol[row] = ZERO ;
02529                 continue ;
02530             }
02531             /*testrow = 0 ; 
02532             while ( lambda > calcol[testrow] )
02533             {
02534                 testrow++ ;
02535             }*/
02536         if (calcolpos[row]==-1) {
02537                if(row >= (*n_rows-1)) calcolpos[row] = calcolpos[row-1];
02538                if(row <  (*n_rows-1)) calcolpos[row] = calcolpos[row+1];
02539         }
02540         testrow = calcolpos[row];
02541 
02542             /*--------------------------------------------------------------------
02543              * at this point calcol[testrow-1] < lambda <= calcol[testrow] 
02544              * now determine the box position in which the polynomial interpolation is carried through.
02545              * the half width of the box is half the number of fit parameters.
02546              * Now we determine the start position of the fitting box and treat
02547              * the special case of being near the edge.
02548              */
02549 
02550             if ( n_params % 2 == 0 )
02551             {
02552                 half_width = (int)(n_params/2) - 1 ;
02553             }
02554             else
02555             {
02556                 half_width = (int)(n_params/2) ;
02557             }
02558 
02559             firstpos   = testrow - half_width ;
02560             if ( firstpos < 0 )
02561             {
02562                 firstpos = 0 ;
02563             }
02564             else if ( firstpos > ((calimage->ly)-n_params) )
02565             {
02566                firstpos = calimage->ly - n_params ;
02567             }
02568             
02569             if ( qfits_isnan(imagecol[firstpos]) )
02570             {
02571                 retimagecol[row] = ZERO ;
02572                 continue ;
02573             }
02574             
02575 
02576             /* we must rescale the x-values (smaller than 1) for the fitting routine */
02577             dif = calcol[firstpos+n_params-1] - calcol[firstpos] ;
02578             for ( i = 0 ; i < n_params ; i++ )
02579             {
02580             x_renorm[i] = (calcol[firstpos + i] - calcol[firstpos]) / dif ;
02581         }
02582         
02583            
02584 
02585             lambda_renorm = ( lambda - calcol[firstpos] ) / dif ; 
02586         
02587         imageptr = &imagecol[firstpos] ;
02588         
02589         flag = 0;
02590         poly=nev_ille(x_renorm, imageptr, n_params-1, lambda_renorm, &flag);
02591         
02592             new_sum += poly ;
02593             retimagecol[row] = poly ; 
02594     }
02595     
02596 
02597         /* now renorm the total flux */
02598         for ( row = 0 ; row < retImage -> ly ; row++ )
02599         {
02600             if ( new_sum == 0. ) new_sum = 1. ;
02601             if ( qfits_isnan(retimagecol[row]) )
02602             {
02603                 retImage->data[col+row*retImage->lx] = ZERO ;
02604             }
02605             else
02606             {
02607                 /* rescaling is commented out because it delivers wrong results
02608                    in case of appearance of blanks or bad pixels */
02609                 retImage->data[col+row*retImage->lx] = retimagecol[row] /* * sum/new_sum*/ ;
02610             }
02611         }
02612 
02613     }
02614     
02615     return retImage ;
02616 }
02617 
02618 
02619 /*___oOo___*/

Generated on Wed Oct 26 13:08:56 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001