focus_det.c

00001 /*----------------------------------------------------------------------------
00002    
00003    File name    :       dark.c
00004    Author       :   A. Modigliani
00005    Created on   :   Sep 17, 2003
00006    Description  :   Focus determination 
00007 
00008  ---------------------------------------------------------------------------*/
00009 
00010 /*----------------------------------------------------------------------------
00011                                 Includes
00012  ---------------------------------------------------------------------------*/
00013 #include "focus_det.h"
00014 
00015 /*----------------------------------------------------------------------------
00016                                 Defines
00017  ---------------------------------------------------------------------------*/
00018 
00019 
00020 /*----------------------------------------------------------------------------
00021                              Function Definitions
00022  ---------------------------------------------------------------------------*/
00023 
00024 /*----------------------------------------------------------------------------
00025    Function     :       focus_det()
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 it does the image reconstruction of a set of offset-corrected, flatfielded, 
00031 bad pixel corrected and slope of the spectra aligned exposures of a point 
00032 source with continuum spectrum. The exposures are taken with different focus 
00033 positions of the detector.
00034 
00035 Afterwards a 2-D Gaussian is fitted to the reconstructed images of the point 
00036 source. The reconstructed images are stored and the sequence of the resulting 
00037 fit parameters (position, amplitude, background, fwhmx, fwhmy, angle of 
00038  major axis) are stored in an ASCII file.
00039 
00040  ---------------------------------------------------------------------------*/
00041 
00042 
00043 int focus_det (cpl_parameterlist* config,cpl_frameset* sof)
00044 {
00045     const char* _id = "focus";
00046     focus_config * cfg ;
00047 
00048     fits_header *  head;
00049 
00050     OneImage * im ;
00051     OneImage * img ;
00052     OneImage * clean_im ;
00053     OneImage * out_image ;
00054     OneCube * cube ;
00055     OneCube * out_cube ;
00056 
00057     int npar =7;
00058     char* name;
00059     char* tbl_name=NULL;
00060 
00061     int i = 0;
00062     int fra=0;
00063     int nslits=0;
00064 
00065     float** slit_edges;
00066     int*  mpar;
00067     double* derv_par;
00068     double* fit_par;
00069     float* correct_dist;
00070     float* vec_x;
00071     float* vec_y;
00072 
00073     double par=0.;
00074     int iters=0;
00075     float* distances;
00076     float fc=0;
00077     char* out_Name_No;
00078     char* gaussplotNameNo;
00079     char* string_par;
00080     char* tmp_string_par;
00081 
00082     FILE * fitpar_file ;
00083     FILE * slit_list ;
00084     FILE * dist_list ;
00085 
00086     int cnt = 0 ;
00087     float edge_x=0;
00088     float edge_y=0;
00089     float tmp_float=0;
00090     cpl_table* tbl_focus=NULL;
00091     cpl_propertylist* tbl_focus_phead=NULL;
00092     cpl_propertylist* tbl_focus_head=NULL;
00093     char* col_name=NULL;
00094     char* tbl_slitpos_name=NULL;
00095     cpl_table* tbl_slitpos=NULL;
00096     char* tbl_distances_name=NULL;
00097     cpl_table* tbl_distances=NULL;
00098     char* tbl_firstcol_name=NULL;
00099     cpl_table* tbl_firstcol=NULL;
00100     int* status=NULL;
00101 
00102     cpl_image* o_img=NULL;
00103     cpl_image* g_img=NULL;
00104     cpl_frameset* stk=NULL;
00105     int check=0;
00106 
00107     /* 
00108        -----------------------------------------------------------------
00109        1) parse the file names and parameters to the dark_config data 
00110           structure cfg 
00111        -----------------------------------------------------------------
00112      */
00113     cfg = parse_cpl_input_focus(config,sof,&stk) ;
00114 
00115     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00116       cpl_msg_error(_id,"error parsing cpl input");
00117       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00118       return -1;
00119     }
00120     if (cfg == NULL)
00121     {
00122         cpl_msg_error (_id,"could not parse cpl input!\n") ;
00123         return -1 ;
00124     }
00125     /* replace the following with some C code
00126 
00127 test = 0
00128 try:
00129     os.environ["CALIB_PATH"]
00130 except:
00131     print "environment variable CALIB_PATH not defined!"
00132     test = 1
00133 
00134 if test == 0:
00135     if poslist[0] == '/':
00136         poslist = os.environ["CALIB_PATH"] + poslist
00137     if firstCol[0] == '/':
00138         firstCol = os.environ["CALIB_PATH"] + firstCol
00139 
00140     */
00141 
00142     if (is_fits_file(cfg->poslist) != 1) {
00143       cpl_msg_error(_id,"poslist not given!");
00144        return -1;
00145     }
00146     out_Name_No = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00147     gaussplotNameNo = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00148 
00149       /* store the found fit parameters in an ASCII file */
00150 
00151       if ( NULL == (fitpar_file = fopen (cfg->fitlist, "w" ) ) )
00152       {
00153         cpl_msg_error(_id,"cannot open %s\n", cfg->fitlist) ;
00154         return -1 ;
00155       }
00156 
00157     col_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00158     tbl_focus = cpl_table_new(cfg->nframes);
00159     for (i=0; i< npar; i++) {
00160       sprintf(col_name,"%s%d","par",i);
00161       cpl_table_new_column(tbl_focus,col_name, CPL_TYPE_DOUBLE);
00162       cpl_table_set_column_format(tbl_focus,col_name,"%14.10f");
00163 
00164     }
00165 
00166     printf("nframes=%d\n",cfg->nframes);
00167     for (fra=0; fra < cfg->nframes; fra++) {
00168 
00169       name = cfg->inFrameList[fra];
00170       if(is_fits_file(name) != 1 ) {
00171         cpl_msg_error(_id,"Input file %s is not FITS",cfg->inFrameList[fra]);
00172         return -1;
00173       }
00174       im = load_image(name);
00175    
00176       if (im == NULL) {
00177          cpl_msg_error(_id, " could not load inFrame\n" );
00178          return -1;
00179       }
00180 
00181 
00182       /* read the fits header to change some standard keywords */
00183       head = fits_read_header(name);
00184 
00185       /*
00186       #------------------------------------------------------------------------
00187       #----------------------RECONSTRUCTION------------------------------------
00188       #------------------------------------------------------------------------
00189       */
00190       clean_im = cleanMeanOfColumns( im, cfg->lo_reject, cfg->hi_reject );
00191       if (clean_im == NULL) {
00192         cpl_msg_error(_id," could not take a clean mean of the spectra!\n");
00193         return -1;
00194       }
00195       destroy_image (im);
00196 
00197       /* select north-south-test or fitting of slitlet edges */
00198       correct_dist = new_float_array(cfg->nslits) ;
00199 
00200       if (cfg->northsouthInd == 0) {
00201 
00202         /*READ TFITS TABLE*/
00203         tbl_slitpos_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00204         strcpy(tbl_slitpos_name,cfg->poslist);
00205         tbl_slitpos = cpl_table_load(tbl_slitpos_name,1,0);
00206 
00207     
00208         for (i =0 ; i< cfg->nslits; i++){
00209       
00210            edge_x=cpl_table_get_float(tbl_slitpos,"start",i,status);
00211            edge_y=cpl_table_get_float(tbl_slitpos,"end",i,status);
00212            array2D_set_value(slit_edges,edge_x,i,0);
00213            array2D_set_value(slit_edges,edge_y,i,1);
00214         }
00215     
00216         cpl_table_delete(tbl_slitpos);
00217         cpl_free(tbl_slitpos_name);
00218 
00219         if(cpl_error_get_code() != CPL_ERROR_NONE) {
00220           cpl_msg_error(_id,"error writing tbl %s",tbl_slitpos_name);
00221           cpl_msg_error(_id,(char* ) cpl_error_get_message());
00222           return -1;
00223     }
00224         slit_edges = new_2Dfloat_array(cfg->nslits, 2) ;
00225 
00226         if ( NULL == (slit_list = fopen (cfg->poslist, "r" ) ) )
00227         {
00228             cpl_msg_error(_id,"cannot open %s\n", cfg->poslist) ;
00229             return -1 ;
00230     }
00231 
00232         cnt = 0 ;
00233         while ( fscanf( slit_list, "%f %f", &vec_x[cnt], &vec_y[cnt] ) != EOF )
00234         {
00235            cnt ++ ;
00236         }
00237 
00238         vec_x=new_float_array(cnt);
00239         vec_y=new_float_array(cnt);
00240 
00241         cnt=0;
00242         while ( fscanf( slit_list, "%f %f", &vec_x[cnt], &vec_y[cnt] ) != EOF )
00243         {
00244            cnt ++ ;
00245         }
00246 
00247 
00248         if (cnt != cfg->nslits) {
00249       cpl_msg_error(_id," wrong file format!\n");
00250       fclose(slit_list);
00251       return -1;
00252     }
00253         fclose(slit_list);
00254 
00255 
00256         /* open the ASCII list of the slitlet positions */
00257         for (i=0; i< nslits; i++) {
00258         array2D_set_value(slit_edges, vec_x[i], i, 0);
00259         array2D_set_value(slit_edges, vec_y[i], i, 1);
00260 
00261     }
00262         fclose(slit_list);
00263 
00264       /* build the data cube */
00265       cube = makeCubeSpi( clean_im, slit_edges, correct_dist );
00266         if (cube == NULL) {
00267       cpl_msg_error(_id," could not construct data cube!\n");
00268       return -1;
00269     }
00270         destroy_2Darray(slit_edges, cfg->nslits);
00271       } else {
00272         distances = new_float_array(cfg->nslits-1) ;
00273 
00274 
00275 
00276         if ( NULL == (dist_list = fopen (cfg->poslist, "r" ) ) )
00277         {
00278             cpl_msg_error(_id,"cannot open %s\n", cfg->poslist) ;
00279             return -1 ;
00280     }
00281   
00282         /*READ TFITS TABLE*/
00283         tbl_distances_name=cfg->poslist;
00284          
00285         tbl_distances = cpl_table_load(tbl_distances_name,1,0);
00286         if(cpl_error_get_code() != CPL_ERROR_NONE) {
00287           cpl_msg_error(_id,"error loadting tbl %s",tbl_distances_name);
00288           cpl_msg_error(_id,(char* ) cpl_error_get_message());
00289           return -1;
00290     }
00291         for (i =0 ; i< cfg->nslits-1; i++){
00292        tmp_float=cpl_table_get_float(tbl_distances,
00293                          "slitlet_distance",i,status);
00294            array_set_value(distances,tmp_float,i);
00295      
00296         }
00297         cpl_table_delete(tbl_distances);
00298         if(cpl_error_get_code() != CPL_ERROR_NONE) {
00299           cpl_msg_error(_id,"error writing tbl %s",tbl_distances_name);
00300           cpl_msg_error(_id,(char* ) cpl_error_get_message());
00301           return -1;
00302     }
00303 
00304     /*
00305         cnt = 0 ;
00306     
00307         while ( fscanf( dist_list, "%f", &vec_x ) != EOF )
00308         {
00309            cnt ++ ;
00310         }
00311         fclose(dist_list);
00312 
00313         vec_x = new_float_array(cnt);
00314         dist_list = fopen (cfg->poslist, "r" );
00315       cnt = 0;
00316         while ( fscanf( dist_list, "%f", &vec_x[cnt] ) != EOF )
00317         {
00318            cnt ++ ;
00319         }
00320 
00321 
00322         if (cnt != (cfg->nslits-1)) {
00323       cpl_msg_error(_id," wrong file format!\n");
00324       fclose(dist_list);
00325       return -1;
00326     }
00327         fclose(dist_list);
00328 
00329     */
00330         /* open the ASCII list of the determined distances */ 
00331     /* 
00332         for (i=0; i< cfg->nslits-1; i++) {
00333       array_set_value(distances, vec_x[i], i);
00334     }
00335     */
00336 
00337         /* build the data cube */
00338        /*READ TFITS TABLE*/
00339         tbl_firstcol_name=cfg->firstCol;
00340         if ((tbl_firstcol=cpl_table_load(tbl_firstcol_name,1,0)) != NULL) {
00341            fc=cpl_table_get_float(tbl_firstcol,"first_col",0,status);
00342        cpl_table_delete(tbl_firstcol);
00343     } else {
00344       cpl_msg_error(_id,"table %s not found! Exit!",tbl_firstcol_name);
00345           return -1 ;
00346     }
00347 
00348    
00349         if(cpl_error_get_code() != CPL_ERROR_NONE) {
00350           cpl_msg_error(_id,"error loading tbl %s",tbl_firstcol_name);
00351           cpl_msg_error(_id,(char* ) cpl_error_get_message());
00352           return -1;
00353     }
00354 
00355     /*
00356         if ( NULL == (fc_list = fopen (cfg->firstCol, "r" ) ) )
00357         {
00358             cpl_msg_error(_id,"cannot open %s\n", cfg->firstCol) ;
00359             return -1 ;
00360     }
00361 
00362         cnt = 0 ;
00363         fscanf( fc_list, "%f", &fc );
00364         fclose(fc_list);
00365     */
00366 
00367         cube = makeCubeDist( clean_im, fc, distances, correct_dist );
00368         if (cube == NULL) {
00369       cpl_msg_error(_id," could not construct a data cube!\n");
00370       return -1;
00371     }
00372         destroy_array (distances);
00373 
00374       }
00375       /*
00376         #----------------------------------------------------------------------
00377         #------------------------FINETUNING------------------------------------
00378         #----------------------------------------------------------------------
00379       */
00380       /* shift the rows of the reconstructed images of the data cube to 
00381            the correct sub pixel position
00382            select the shift method: polynomial interpolation, FFT or cubic 
00383            spline interpolation
00384       */
00385       if (strcmp(cfg->method,"P") == 0) {
00386       out_cube = fineTuneCube( cube, correct_dist, cfg->order );
00387       if (out_cube == NULL) {
00388             cpl_msg_error(_id," could not fine tune the data cube\n");
00389             return -1;
00390       }
00391       } else if (strcmp(cfg->method, "F") == 0) {
00392         for (i=0; i< cfg->nslits; i++) {
00393           array_set_value(correct_dist, -array_get_value(correct_dist,i), 
00394                               i);
00395         }
00396         out_cube = fineTuneCubeByFFT( cube, correct_dist );
00397         if (out_cube == NULL) {
00398           cpl_msg_error(_id," could not fine tune the data cube\n");
00399           return -1;
00400         }
00401       } else if ( strcmp(cfg->method, "S") == 0) {
00402       out_cube = fineTuneCubeBySpline( cube, correct_dist );
00403       if (out_cube == NULL) {
00404             cpl_msg_error(_id," could not fine tune the data cube\n");
00405             return -1;
00406       }
00407       } else {
00408       cpl_msg_error(_id," wrong method indicator given!");
00409       return -1;
00410       }
00411       /* convert cube format to image format */
00412       out_image = medianCube(out_cube);
00413       if (out_image == NULL) {
00414       cpl_msg_error(_id," could not do medianCube()\n");
00415       return -1;
00416       }
00417 
00418       if (strstr(cfg->outName, "." ) != NULL ) {
00419       sprintf(out_Name_No, "%s%d%s", 
00420           get_rootname(cfg->outName), fra,strstr(cfg->outName,"."));
00421       } else {
00422       sprintf(out_Name_No, "%s%d", cfg->outName, fra);
00423       }
00424 
00425       o_img=cpl_image_wrap_float(out_image->lx,out_image->ly,out_image->data);
00426     
00427 
00428 
00429       if(-1 == sinfoni_pro_dump_ima(o_img,stk,sof,out_Name_No,
00430                 PRO_FOCUS,_id,config)) {
00431                 cpl_msg_error(_id,"cannot dump ima %s", out_Name_No);
00432                 destroy_array(correct_dist);
00433                 destroy_image(clean_im);
00434                 destroy_cube(cube);
00435                 destroy_cube (out_cube);
00436                 destroy_image (out_image);
00437                 fits_header_destroy(head);
00438                 cpl_free(col_name);
00439                 fclose(fitpar_file);
00440                 cpl_free(out_Name_No);
00441                 cpl_free(gaussplotNameNo);
00442                 free_focus(cfg);
00443         return -1;
00444       }
00445  
00446       /*
00447         #----------------------------------------------------------------
00448         #------------------------------GAUSS2DFIT------------------------
00449         #----------------------------------------------------------------
00450        */
00451        /* allocate memory and associate the mask parameters with 
00452              given values */
00453 
00454       fit_par  = new_double_array(npar);
00455       derv_par = new_double_array(npar);
00456       mpar     = new_int_array(npar);
00457 
00458 
00459       intarray_set_value(mpar, cfg->mpar0, 0);
00460       intarray_set_value(mpar, cfg->mpar1, 1);
00461       intarray_set_value(mpar, cfg->mpar2, 2);
00462       intarray_set_value(mpar, cfg->mpar3, 3);
00463       intarray_set_value(mpar, cfg->mpar4, 4);
00464       intarray_set_value(mpar, cfg->mpar5, 5);
00465       intarray_set_value(mpar, cfg->mpar6, 6);
00466         
00467       iters = fit2dGaussian( out_image, fit_par, derv_par, 
00468                            mpar, cfg->llx, cfg->lly, 
00469                            cfg->halfbox_x, cfg->halfbox_y,&check );
00470       if (iters < 0) {
00471         cpl_msg_error(_id," could not fit 2d-Gaussian\n");
00472         return -1;
00473       }
00474 
00475       /* here there was open file */
00476 
00477       string_par = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00478       tmp_string_par = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00479 
00480       sprintf(string_par,"%s","");
00481       for (i=0; i< npar; i++) {
00482         sprintf(col_name,"%s%d","par",i);
00483         par = doublearray_get_value(fit_par, i);
00484         sprintf(tmp_string_par,"%14.10f%s",par," ");
00485        
00486         strcat(string_par,tmp_string_par);
00487         cpl_table_set_double(tbl_focus,col_name,fra,par);
00488       }
00489       fprintf(fitpar_file, "%s \n", string_par);
00490       cpl_free(string_par);
00491       cpl_free(tmp_string_par);
00492 
00493       if (cfg->plotGaussInd == 1) {
00494         img = plotGaussian( out_image, fit_par ); 
00495         if (img == NULL) {
00496       cpl_msg_error(_id," could not do plotGaussian()\n");
00497       return -1;
00498     }
00499 
00500         if (strstr(cfg->gaussplotName, "." ) != NULL ) {
00501       sprintf(gaussplotNameNo, "%s%d%s", 
00502           get_rootname(cfg->gaussplotName), fra,
00503           strstr(cfg->gaussplotName,"."));
00504 
00505         } else {
00506       sprintf(gaussplotNameNo, "%s%d", cfg->gaussplotName, fra);
00507         }
00508 
00509       g_img=cpl_image_wrap_float(img->lx,img->ly,img->data);
00510 
00511 
00512       if(-1 == sinfoni_pro_dump_ima(g_img,stk,sof,gaussplotNameNo,
00513                 PRO_FOCUS_GAUSS,_id,config)) {
00514                 cpl_msg_error(_id,"cannot dump ima %s", gaussplotNameNo);
00515 
00516                 destroy_intarray(mpar);
00517                 destroy_doublearray(fit_par);
00518                 destroy_doublearray(derv_par);
00519                 destroy_image(img);
00520 
00521                 destroy_array(correct_dist);
00522                 destroy_image(clean_im);
00523                 destroy_cube(cube);
00524                 destroy_cube (out_cube);
00525                 destroy_image (out_image);
00526                 fits_header_destroy(head);
00527                 cpl_free(col_name);
00528                 fclose(fitpar_file);
00529                 cpl_free(out_Name_No);
00530                 cpl_free(gaussplotNameNo);
00531                 free_focus(cfg);
00532         return -1;
00533       }
00534     
00535 
00536         destroy_image(img);
00537     
00538       }
00539 
00540       /* ---free memory---- */
00541       destroy_intarray(mpar);
00542       destroy_doublearray(fit_par);
00543       destroy_doublearray(derv_par);
00544       destroy_array(correct_dist);
00545       destroy_image(clean_im);
00546       destroy_cube(cube);
00547       destroy_cube (out_cube);
00548       destroy_image (out_image);
00549       fits_header_destroy(head);
00550 
00551     }   
00552     tbl_name=cfg->fitlist;
00553     cpl_table_save(tbl_focus,tbl_focus_phead,tbl_focus_head,tbl_name,0);
00554     cpl_table_delete(tbl_focus);
00555     cpl_free(col_name);
00556     fclose(fitpar_file);
00557     cpl_free(out_Name_No);
00558     cpl_free(gaussplotNameNo);
00559     free_focus(cfg);
00560 
00561 
00562     return 0 ;
00563 }
00564 

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