wave_cal.c

00001 /*----------------------------------------------------------------------------
00002  
00003      File name    :     wave_cal.c
00004    Author       :   A. Modigliani
00005    Created on   :   Sep 17, 2003
00006    Description  : 
00007 
00008  ---------------------------------------------------------------------------*/
00009 
00010 /*----------------------------------------------------------------------------
00011                                 Includes
00012  ---------------------------------------------------------------------------*/
00013 #include "wave_cal.h"
00014 
00015 /*----------------------------------------------------------------------------
00016                                 Defines
00017  ---------------------------------------------------------------------------*/
00018 
00019 
00020 /*----------------------------------------------------------------------------
00021                              Function Definitions
00022  ---------------------------------------------------------------------------*/
00023 
00024 /*----------------------------------------------------------------------------
00025    Function     :       wave_cal()
00026    In           :       ini_file: file name of according .ini file
00027    Out          :       integer (0 if it worked, -1 if it doesn't) 
00028    Job          :
00029 
00030 
00031   Normal method:
00032 
00033   does the wavelength calibration and the fitting of the slitlet edge 
00034   positions (ASCII file 32 x 2 values) if wished
00035   produces an array of the bcoefs and of the fit parameters if wished and a 
00036   wavelength calibration map input is an emission line frame and a line list
00037 
00038 
00039   o searching for lines by cross correlation with a line list
00040   o Gaussian fitting of emission lines in each column->positions of the lines->
00041     resulting fit parameters can be stored in an ASCII file
00042   o Fitting of a polynomial to the line positions for each column
00043   o Smoothing: fitting of each polynomial coefficient by another polynomial
00044     across the whole frame -> resulting polynomial coefficients can be stored 
00045     in an ASCII file.
00046   o Wavelength calibration map (micron value for each frame pixel) can be
00047     produced by using these coefficients and a cross correlation to the
00048     original frame
00049 
00050   o The slitlet edge positions can be fitted:
00051     1) Automatically (not really stable) or by using guess edge positions
00052     2) By using a Boltzmann or a linear slope function
00053 
00054   Slit method:
00055 
00056   does the wavelength calibration and the fitting of the slitlet edge 
00057   positions (ASCII file 32 x 2 values) if wished produces a list of the fit 
00058   parameters and of the smoothed coefficients if wished and a wavelength 
00059   calibration map input is an emission line frame and a line list
00060 
00061   o Does the same as other method but smoothes the found polynomial 
00062     coefficients within each slitlet and not over the whole frame.
00063 
00064   o Produces always a wavelength calibration map and does not crosscorrelate.
00065 
00066  ---------------------------------------------------------------------------*/
00067 
00068 
00069 int wave_cal (cpl_parameterlist* config, cpl_frameset* sof) {
00070     
00071     const char* _id = "wave_cal";
00072 
00073     wave_config * cfg ;
00074     OneImage *  map ;
00075 
00076     OneImage * im ;
00077     FitParams** par;
00078     float* wave;
00079     float* intens;
00080     float** bcoefs;
00081 
00082     int check = 0;
00083     int lx = 0;
00084     int ly = 0;
00085     int n_lines=0;
00086     int i = 0;
00087     int j = 0;
00088     int n = 0;
00089     int* n_found_lines;
00090     int* sum_pointer;
00091     int** row_clean;
00092     float** wavelength_clean;
00093 
00094     int sum=0;
00095     int er=0;
00096     float a=0.;
00097 
00098     float val1=0.;
00099     float val2=0.;
00100     float val3=0.;
00101     float val4=0.;
00102     float val5=0.;
00103 
00104     float** val_x;
00105     float** val_y;
00106 
00107     float* vec_x;
00108     float* vec_y;
00109     int cnt = 0 ;
00110 
00111     float** slit_pos;
00112     char* tmp_str;
00113     int readsum=0;
00114     int fit=0;
00115     float** val;
00116     FILE * fp ;
00117     FILE * bcoefs_file;
00118     FILE * bcoefs_list;
00119     FILE * par_list;
00120     FILE * slit_list;
00121     FILE * slitpos_file;
00122     cpl_frameset* raw=NULL;
00123 
00124     /* 
00125        -----------------------------------------------------------------
00126        1) parse the file names and parameters to the ns_config data 
00127           structure cfg
00128        -----------------------------------------------------------------
00129      */
00130 
00131     cfg = parse_cpl_input_wave(config,sof,&raw) ;
00132     if (cfg == NULL)
00133     {
00134         cpl_msg_error ("wave_cal:","Error could not parse cpl input!") ;
00135         return -1 ;
00136     }
00137     /* Convert this python code into C
00138 
00139 test = 0
00140 try:
00141     os.environ["CALIB_PATH"]
00142 except:
00143     print "environment variable CALIB_PATH not defined!"
00144     test = 1
00145 
00146 if test == 0:
00147     if lineList[0] == '/':
00148         lineList = os.environ["CALIB_PATH"] + lineList
00149     if slitposGuessName[0] == '/':
00150         slitposGuessName = os.environ["CALIB_PATH"] + slitposGuessName
00151     if calibIndicator == 0 and wavemapInd == 1:
00152         if coeffsName[0] == '/':
00153             coeffsName = os.environ["CALIB_PATH"] + coeffsName
00154     if slitposIndicator == 1:
00155         if calibIndicator != 1 and estimateIndicator != 1:
00156             if paramsList[0] == '/':
00157                 paramsList = os.environ["CALIB_PATH"] + paramsList
00158 
00159 
00160 
00161     */
00162 
00163     if(is_fits_file (cfg->inFrame) != 1) {
00164       cpl_msg_error(_id,"Input file %s is not FITS",cfg->inFrame);
00165       return -1;
00166     }
00167     
00168     if (cfg->slitposIndicator == 1 && cfg->estimateIndicator == 1) {
00169       if (is_fits_file(cfg->slitposGuessName) != 1) {
00170         cpl_msg_error("wave_cal:","slitlet position guess list not given!");
00171         exit(-1);
00172       }
00173     }
00174 
00175     if (cfg->calibIndicator == 0 && cfg->wavemapInd == 1) {
00176       if (is_fits_file(cfg->coeffsName) != 1) {
00177         cpl_msg_error("wave_cal:","coefficients list not given!");
00178         exit(-1);
00179 
00180       }
00181     }
00182     if (cfg->slitposIndicator == 1) {
00183       if (cfg->calibIndicator != 1 && cfg->estimateIndicator != 1) {
00184         if (is_fits_file(cfg->paramsList) != 1) {
00185        cpl_msg_error("wave_cal:","parameter list not given!");
00186            exit(-1);
00187     }
00188       }
00189     }
00190 
00191     /* ---load the emission line frame--- */
00192     im = load_image(cfg->inFrame);
00193     if (im == NULL) {
00194       cpl_msg_error ("wave_cal","could not load image");
00195       exit(-1);
00196     }
00197 
00198     lx = im->lx;
00199     ly = im->ly;
00200 
00201     if (cfg->calibIndicator == 1 || cfg->wavemapInd == 1) {
00202       if (is_fits_file(cfg->lineList) != 1) {
00203         cpl_msg_error("wave_cal:","linelist not given!");
00204         exit(-1);
00205       }
00206     }
00207 
00208     /* ---open the line list and read the number of lines--- */
00209     /* The following to be adjusted from py2C conversion */
00210     if ( NULL == (fp = fopen (cfg->lineList, "r" ) ) )
00211     {
00212         cpl_msg_error("wave_cal:","cannot open %s", cfg->lineList) ;
00213         return -1 ;
00214     }
00215 
00216     cnt = 0 ;
00217     while ( fscanf( fp, "%f %f", &vec_x[cnt], &vec_y[cnt] ) != EOF )
00218     {
00219         cnt ++ ;
00220     }
00221     n= cnt ;
00222 
00223     vec_x = new_float_array (n);
00224     vec_y = new_float_array (n);
00225 
00226     cnt=0;
00227     while ( fscanf( fp, "%f %f", &vec_x[cnt], &vec_y[cnt] ) != EOF )
00228     {}
00229 
00230 
00231     /* #---allocate memory for the line list--- */
00232     wave   = new_float_array (n);
00233     intens = new_float_array (n);
00234     if (wave == NULL || intens == NULL) {
00235       cpl_msg_error ( "wave_cal","could not allocate memory for the line list values" );
00236       exit(-1);
00237     }
00238 
00239     /* #---read the line list--- */
00240     n_lines = readList(cfg->lineList, wave, intens);
00241     if (-1 == n_lines) { 
00242        cpl_msg_error("wave_cal","readList failed!");
00243        exit(-1);
00244     }
00245 
00246     /* this is probably to be converted from python to C */
00247     fclose(fp);   
00248     bcoefs  = new_2Dfloat_array(cfg->nrDispCoefficients, cfg->nrCoefCoefficients);
00249 
00250     /*----------------------------------------------------------------
00251      *---------------------FINDLINES----------------------------------
00252      *----------------------------------------------------------------
00253      */
00254     /*
00255       if not yet done: do the wavelength calibration, that means: 
00256       find the dispersion relation and parameterize its coefficients
00257      */
00258 
00259     if (cfg->calibIndicator == 1) {
00260       /* allocate memory */
00261       n_found_lines    = new_int_array(lx); 
00262       row_clean        = new_2Dint_array(lx, n_lines);
00263       wavelength_clean = new_2Dfloat_array(lx, n_lines);
00264       sum_pointer      = new_int_array(1);
00265  
00266       /* find the emission lines in each image column */
00267       intarray_set_value(sum_pointer, 0, 0);
00268       check = findLines(im, wave, intens, n_lines, row_clean, 
00269                         wavelength_clean, cfg->guessBeginWavelength, 
00270                         cfg->guessDispersion1, cfg->guessDispersion2,
00271                         cfg->mindiff, cfg->halfWidth,
00272                         n_found_lines, cfg->sigma, sum_pointer );
00273       if (-1 == check) { 
00274         cpl_msg_error ("wave_cal","findLines failed!");
00275         exit(-1);
00276       }
00277 
00278 
00279       /*---------------------------------------------------------------------
00280        *-----------------------WAVE_CALIB-------------------------------------
00281        *---------------------------------------------------------------------
00282        */
00283 
00284       sum = intarray_get_value(sum_pointer,0);
00285 
00286       /* allocate memory for the fit parameters */
00287       par = newFitParams( sum );
00288       if (par == NULL) {
00289         cpl_msg_error("wave_cal","newFitParams failed!");
00290         exit(-1);
00291       }
00292 
00293       /* 
00294          fit each line, make a polynomial fit and fit the resulting fit 
00295          coefficients across the columns
00296        */
00297 
00298       er = wavelengthCalibration(im, par, bcoefs, wave, n_lines, row_clean, 
00299                                wavelength_clean, n_found_lines, 
00300                    cfg->guessDispersion1, 
00301                                cfg->halfWidth, cfg->minAmplitude, 
00302                                cfg->maxResidual, cfg->fwhm, 
00303                                cfg->nrDispCoefficients, 
00304                                cfg->nrCoefCoefficients, 
00305                                cfg->sigmaFactor, 
00306                                cfg->pixel_tolerance);
00307       if (-1 == er) { 
00308          cpl_msg_error ("wave_cal","wavelengthCalibration failed!");
00309          exit(-1);
00310       }
00311 
00312       /* store the resulting polynomial fit coefficients in an ASCII 
00313          file if wished
00314        */
00315 
00316       tmp_str = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char));
00317       if (cfg->writeCoeffsInd == 1) {
00318 
00319     if ( NULL == (fp = fopen ( cfg->coeffsName, "w" ) ) ) 
00320          {
00321            cpl_msg_error("wave_cal:","cannot open %s", cfg->coeffsName) ;
00322            return -1 ;
00323          }
00324          for (i = 0; i< cfg->nrDispCoefficients; i++) {
00325         for (j=0; j< cfg->nrCoefCoefficients; j++) {
00326            a = array2D_get_value(bcoefs, i, j);
00327                sprintf(tmp_str,"%f",a);
00328                strcat(tmp_str," \n");
00329               fprintf(bcoefs_file,"%s",tmp_str);
00330         }
00331 
00332             if ( 0 == fclose ( bcoefs_file ) )
00333             {
00334                cpl_msg_error("wave_cal:","cannot close %s", cfg->coeffsName) ;
00335                return -1 ;
00336         }
00337      }
00338          /* this is python...replaced by C lines above
00339          bcoefs_file = open(cfg->coeffsName, "w");  
00340          for (i = 0; i< cfg->nrDispCoefficients; i++) {
00341         for (j=0; j< cfg->nrCoefCoefficients; j++) {
00342            a = array2D_get_value(bcoefs, i, j);
00343            b.append(str(a)+" ");
00344            bcoefs_file.writelines(b);
00345            bcoefs_file.write("\n");
00346         }
00347             bcoefs_file.close();
00348          }
00349          */
00350       }
00351 
00352       /* store the resulting Gaussian fit parameters in an ASCII 
00353          file if wished
00354        */
00355       if (cfg->writeParInd == 1) {
00356           dumpFitParamsToAscii(par, cfg->paramsList);
00357       }
00358       /* # free memory */
00359       destroy_2Darray ( wavelength_clean, lx );
00360       destroy_2Dintarray (row_clean, lx);
00361       destroy_intarray ( n_found_lines );
00362       destroy_intarray ( sum_pointer );
00363 
00364     }
00365       /*----------------------------------------------------------------------
00366        *-------------------WAVEMAP--------------------------------------------
00367        *----------------------------------------------------------------------
00368        */
00369       /* now do the cross correlation and produce a wavelength map--- */
00370 
00371 
00372     if (cfg->wavemapInd == 1) {
00373       /* #read the parameterized dispersion relation */
00374 
00375 
00376        if ( NULL == (bcoefs_list = fopen (cfg->coeffsName, "r" ) ) )
00377        {
00378            cpl_msg_error("wave_cal:","cannot open %s", cfg->coeffsName) ;
00379            return -1 ;
00380        }
00381 
00382        for (i = 0; i< cfg->nrDispCoefficients; i++) {
00383          for (j = 0; j< cfg->nrCoefCoefficients; j++) {
00384          fscanf( bcoefs_list, "%f", &val[i][j]);
00385          array2D_set_value (bcoefs, val[i][j], i, j);
00386      }
00387        }
00388        fclose(bcoefs_list);
00389 
00390        map = waveMap ( im, bcoefs, cfg->nrDispCoefficients, 
00391                                    cfg->nrCoefCoefficients, 
00392                     wave, intens, n_lines, cfg->magFactor );
00393       if ( map == NULL ) {
00394         cpl_msg_error("wave_cal:","waveMap failed!");
00395         exit(-1);
00396       }
00397       save_image_to_fits( map, cfg->outName, -32 );
00398 
00399       /* # ---free memory--- */
00400       destroy_image(map);
00401 
00402     }
00403       /*-------------------------------------------------------------------
00404        *-------------------SLITFITS----------------------------------------
00405        *-------------------------------------------------------------------
00406        */
00407       /* #--fit the slitlet edge positions if desired-- */
00408 
00409       if (cfg->slitposIndicator == 1) {
00410     if (cfg->calibIndicator != 1 && cfg->estimateIndicator != 1) {
00411       /* read the first integer to determine the number of fit 
00412              parameters to allocate
00413       */
00414 
00415 
00416           if ( NULL == (par_list = fopen (cfg->paramsList, "r" ) ) )
00417           {
00418               cpl_msg_error("wave_cal:","cannot open %s", cfg->paramsList) ;
00419               return -1 ;
00420           }
00421 
00422           cnt = 0 ;
00423           while ( fscanf( par_list, "%f %f %f %f %f", 
00424             &val1,&val2,&val3,&val4,&val5 ) != EOF )
00425           {
00426               cnt ++ ;
00427           }
00428           n= cnt ;
00429       readsum = val5;
00430       fclose(par_list);
00431 
00432 
00433 
00434 
00435 
00436 
00437           /* allocate memory for the fit parameters and read them 
00438              from an ASCII file
00439       */
00440       par = newFitParams( readsum );
00441       if (par == NULL ) {
00442             cpl_msg_error("wave_cal","newFitParams failed!");
00443             exit(-1);
00444       }
00445       dumpAsciiToFitParams(par, cfg->paramsList);
00446 
00447     }
00448     /* allocate memory for the slitlet position array */
00449     slit_pos = new_2Dfloat_array(32,2);
00450     /* now decide which fit model function you want to use by 
00451              read a given character
00452      */
00453     if (cfg->fitBoltzIndicator == 1) {
00454         if (cfg->estimateIndicator == 1) {
00455           /* open the ASCII list of the slitlet positions----- */
00456 
00457 
00458               if ( NULL == (slit_list = fopen (cfg->slitposGuessName, "r" ) ) )
00459               {
00460                   cpl_msg_error("wave_cal:","cannot open %s", cfg->slitposGuessName) ;
00461                   return -1 ;
00462               }
00463 
00464             
00465           for (i =0; i< 32; i++){
00466                    fscanf(slit_list, "%f %f", 
00467               &val_x[i][0],&val_y[i][1] );
00468            
00469           array2D_set_value(slit_pos, val_x[i][0], i, 0);
00470           array2D_set_value(slit_pos, val_y[i][1], i, 1);
00471 
00472           }
00473           fclose(slit_list);
00474          
00475 
00476               fit = fitSlitsBoltzWithEstimate( im, slit_pos, cfg->boxLength, 
00477                          cfg->yBox, cfg->diffTol, 
00478                                              cfg->loPos, cfg->hiPos );
00479               if (fit < 0){
00480             cpl_msg_error ( "wave_cal:","fitSlitsBoltz failed" );
00481             exit(-1);
00482           }
00483 
00484         } else {
00485 
00486           fit = fitSlitsBoltz( im, par, slit_pos, 
00487                                    cfg->boxLength, cfg->yBox, 
00488                                    cfg->diffTol );
00489 
00490         }
00491     } else if ( cfg->fitEdgeIndicator == 1) {
00492           if (cfg->estimateIndicator == 1) {
00493         /* open the ASCII list of the slitlet positions-- */
00494 
00495                 if (NULL == (slit_list = fopen (cfg->slitposGuessName, "r" )))
00496                 {
00497                   cpl_msg_error("wave_cal:","cannot open %s", cfg->slitposGuessName) ;
00498                   return -1 ;
00499                 }
00500 
00501             
00502             for (i =0; i< 32; i++){
00503                    fscanf(slit_list, "%f %f", 
00504               &val_x[i][0],&val_y[i][1] );
00505            
00506           array2D_set_value(slit_pos, val_x[i][0], i, 0);
00507           array2D_set_value(slit_pos, val_y[i][1], i, 1);
00508 
00509             }
00510             fclose(slit_list);
00511 
00512 
00513                 fit = fitSlitsEdgeWithEstimate( im, slit_pos, cfg->boxLength, 
00514                         cfg->yBox, cfg->diffTol, 
00515                                                 cfg->loPos, cfg->hiPos );
00516         if (fit < 0) {
00517           cpl_msg_error ( "wave_cal:","fitSlitsBoltz failed" );
00518           exit(-1);
00519         }
00520           } else {
00521                   fit = fitSlitsEdge( im, par, slit_pos, 
00522                                       cfg->boxLength, cfg->yBox, 
00523                                       cfg->diffTol );
00524           if (fit < 0) {
00525             cpl_msg_error ( "wave_cal:","fitSlitsEdge failed" );
00526             exit(-1);
00527           } 
00528 
00529           }
00530     } else {
00531           cpl_msg_error ( "wave_cal:","no indication of desired fit function given, will exit" );
00532       exit(-1);
00533 
00534     }
00535         
00536     /* store the resulting sitlet positions in an ASCII file */
00537 
00538         if (NULL == (slitpos_file = fopen (cfg->slitposName, "w" )))
00539         {
00540             cpl_msg_error("wave_cal:","cannot open %s", cfg->slitposName) ;
00541             return -1 ;
00542         }
00543 
00544             
00545     for (i =0; i< 32; i++){
00546             fprintf(slitpos_file, "%f %s %f\n", 
00547                    &slit_pos[i][0]," ", &slit_pos[i][1] );
00548     }
00549     fclose(slitpos_file);
00550     /* free memory */
00551     destroy_2Darray ( slit_pos, 32 );
00552       }
00553     
00554       
00555    
00556     
00557 
00558         /* #-----free the rest memory---------------------- */
00559 if ((cfg->slitposIndicator == 1 &&  cfg->estimateIndicator != 1) || 
00560      cfg->calibIndicator == 1){
00561   destroyFitParams(par);
00562 }
00563  if (cfg->calibIndicator == 1 || cfg->wavemapInd == 1) {
00564    destroy_array(intens);
00565    destroy_array(wave);
00566    destroy_2Darray ( bcoefs, cfg->nrDispCoefficients );
00567  }
00568  destroy_image ( im );
00569  wavecal_free(cfg);
00570  return 0;
00571 
00572 }

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