nst.c

00001 /*----------------------------------------------------------------------------
00002  
00003      File name    :     spiffi_north_south_test.c
00004    Author       :   A. Modigliani
00005    Created on   :   Sep 17, 2003
00006    Description  : 
00007 
00008  * nsh.h
00009  * Result of a north-south test exposure are 32 continuum spectra of a 
00010  * pinhole that means one spectrum in each slitlet at the same spatial 
00011  * position.
00012  * Each spectrum is fitted in sp[atial direction by a Gaussian to get the 
00013  * sub-pixel positions for each row.
00014  *
00015  * Then the distances are determined in each row and averaged
00016  *
00017  * Result: are distances of each slitlet from each other => 31 values stored 
00018  *  in an ASCII file this Python script needs a frame of a pinhole source with 
00019  *  a continuous spectrum, that is shifted exactly perpendicular to the 
00020  *  slitlets. It fits the spectra in spatial direction by a Gaussian fit 
00021  *  function and therefore determines the sub-pixel position of the source.
00022  *
00023  *  Then the distances of the slitlets from each other are determined and 
00024  *   saved in an ASCII list. 
00025  *
00026 
00027 
00028  ---------------------------------------------------------------------------*/
00029 
00030 /*----------------------------------------------------------------------------
00031                                 Includes
00032  ---------------------------------------------------------------------------*/
00033 #include "nst.h"
00034 #include "sinfoni_pro_save.h"
00035 #include "ns_ini_by_cpl.h"
00036 
00037 /*----------------------------------------------------------------------------
00038                                 Defines
00039  ---------------------------------------------------------------------------*/
00040 
00041 
00042 /*----------------------------------------------------------------------------
00043                              Function Definitions
00044  ---------------------------------------------------------------------------*/
00045 
00046 /*----------------------------------------------------------------------------
00047    Function     :       north_south_test()
00048    In           :       ini_file: file name of according .ini file
00049    Out          :       integer (0 if it worked, -1 if it doesn't) 
00050    Job          : Result of a north-south test exposure are 32 continuum 
00051                   spectra of a pinhole  that means one spectrum in each 
00052                   slitlet at the same spatial position.
00053 
00054                   Each spectrum is fitted in sp[atial direction by a Gaussian
00055                   to get the sub-pixel positions for each row.
00056 
00057                   Then the distances are determined in each row and averaged
00058 
00059                   Result: are distances of each slitlet from each other =>
00060                           31 values stored in an ASCII file
00061 
00062 
00063  ---------------------------------------------------------------------------*/
00064 int nst (const char* plugin_id, cpl_parameterlist* config, cpl_frameset* sof)
00065 {
00066     const char* _id = "nst";
00067     ns_config * cfg ;
00068     OneImage ** list_images ;
00069     OneImage ** list_object ;
00070     OneImage ** list_off;
00071 
00072     OneImage * im_on ;
00073     OneImage * im_on_sub ;
00074     OneImage * im_on_ind ;
00075     OneImage * im_on_gauss ;
00076     OneImage * im_mask ;
00077     OneImage * im_off ;
00078     OneCube  * cube_object ;
00079     OneCube  * cube_off ;
00080     cpl_image* im_on_ind_imgw=NULL;
00081     char* name;
00082     char* tbl_name=NULL;
00083     int typ;
00084     int i = 0;
00085 
00086     int nob =0;
00087     int nof =0;
00088 
00089     float* distances=NULL;
00090     cpl_image* im_on_ind_img=NULL ;
00091     cpl_table* tbl_dist=NULL;
00092 
00093     cpl_vector* qc_dist=NULL;
00094 
00095     double  qc_dist_mean=0;
00096     double  qc_dist_stdev=0;
00097     cpl_frameset* raw=NULL;
00098 
00099     cpl_table* qclog_tbl=NULL;
00100     char* key_value=NULL;
00101     char* key_name=NULL;
00102 
00103     /* 
00104        -----------------------------------------------------------------
00105        1) parse the file names and parameters to the ns_config data 
00106           structure cfg
00107        -----------------------------------------------------------------
00108      */
00109 
00110     /* 
00111       parse the file names and parameters to the ns_config data structure cfg 
00112      */
00113     cpl_msg_info(_id,"Parse cpl input");
00114      raw=cpl_frameset_new();
00115     cfg = parse_cpl_input_ns(config,sof,&raw) ;
00116     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00117       cpl_msg_error(_id,"error parsing cpl input");
00118       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00119       cpl_frameset_delete(raw);
00120       ns_free (cfg);
00121       return -1;
00122     }
00123 
00124     if (cfg == NULL)
00125     {
00126         cpl_msg_error (_id,"could not parse cpl input!\n") ;
00127         cpl_frameset_delete(raw);
00128         return -1 ;
00129 
00130     }
00131 
00132 
00133     if (cfg->maskInd == 1) {
00134       if(is_fits_file(cfg->mask) != 1) {
00135     cpl_msg_error(_id,"Input file %s is not FITS",cfg->mask); 
00136     ns_free (cfg);
00137     cpl_frameset_delete(raw);
00138         return -1;
00139       }
00140     }
00141     /*
00142      --------------------------------------------------------------------
00143      stack the frames in data cubes and take the clean mean of all frames
00144      --------------------------------------------------------------------
00145     */
00146     /* allocate memory for lists of object and off-frames */
00147     cpl_msg_info(_id,"stack the frames in data cubes");
00148     list_object = new_image_list(cfg->nobj);
00149       if (list_object == NULL) {
00150         cpl_msg_error(_id,"could not allocate memory\n");
00151     ns_free (cfg);
00152     cpl_frameset_delete(raw);
00153         return -1;
00154       }
00155 
00156     if (cfg->noff > 0 ) {
00157       list_off = new_image_list(cfg->noff);
00158       if (list_off == NULL) {
00159         cpl_msg_error(_id,"could not allocate memory\n");
00160         destroy_list(list_object);
00161     ns_free (cfg);
00162     cpl_frameset_delete(raw);
00163         return -1;
00164       }
00165     }
00166 
00167 
00168  
00169 /*
00170       #build different image lists for the different cases----
00171 */
00172     cpl_msg_info(_id,"build different image lists for the different cases");
00173     nob = 0;
00174     nof = 0;
00175     list_images = (OneImage**) cpl_calloc(cfg -> nframes, sizeof(OneImage*));
00176 
00177     for ( i = 0; i< cfg->nframes; i++) {
00178        name = cfg->framelist[i];
00179       if(is_fits_file(name) != 1) {
00180         cpl_msg_error(_id,"Input file %s is not FITS",name);
00181         destroy_list(list_object);
00182     ns_free (cfg);
00183     cpl_frameset_delete(raw);
00184         return -1;
00185       }
00186        list_images[i]=load_image( name );
00187 /*
00188        #---read the fits header to change the gain and noise parameter-----
00189 */
00190     }
00191 
00192 
00193 for (i=0; i< cfg->nframes; i++){
00194   typ = intarray_get_value( cfg->frametype, i );
00195   if (list_images[i] == NULL) {
00196     cpl_msg_error(_id,"could not load image\n");
00197     return -1;
00198   }
00199   if (typ == 1) {
00200       insert_image( list_object, list_images[i], nob );
00201       nob = nob + 1;
00202   } else {
00203       insert_image( list_off, list_images[i], nof);
00204       nof = nof + 1;
00205   }
00206 }
00207 
00208 
00209 
00210 
00211 if (cfg->noff != nof || cfg->nobj != nob ){       
00212       cpl_msg_error(_id,"something wrong with the number of the different types of frames");
00213       return -1;
00214 }
00215 
00216 
00217  cube_object = list_make_cube(list_object, nob);
00218  if (cube_object == NULL){
00219       cpl_msg_error (_id,"could not create data cube!");
00220       return -1;
00221  }
00222  destroy_list(list_object);
00223 
00224  if (nof > 0) {
00225    cube_off = list_make_cube(list_off, nof);
00226    if (cube_off == NULL) {
00227         cpl_msg_error (_id,"could not create data cube!");
00228         return -1;
00229    }
00230    destroy_list(list_off);
00231  }
00232 
00233 
00234 
00235 
00236   
00237 /*
00238 #---take the average of the different cubes -------------
00239 */
00240  cpl_msg_info(_id,"take the average of the different cubes");
00241  im_on = average_with_rejection( cube_object, 
00242                                     cfg->loReject, cfg->hiReject );
00243  if (im_on == NULL ) {
00244       cpl_msg_error (_id, "average_with_rejection failed\n" );
00245       return -1;
00246  }
00247  if (cfg->noff != 0) {
00248       im_off = average_with_rejection( cube_off, 
00249                                        cfg->loReject, cfg->hiReject );
00250       if (im_off == NULL) {
00251         cpl_msg_error(_id,"average_with_rejection failed\n" );
00252         return -1;
00253       }
00254       destroy_cube(cube_off);
00255  }
00256  destroy_cube(cube_object);
00257 
00258 
00259 
00260 /*
00261 #finally, subtract off from on frames and store the result in the object cube
00262 */
00263 
00264  if (cfg->noff != 0) {
00265    cpl_msg_info(_id,"subtract off from on frames");
00266       im_on_sub = sub_image(im_on, im_off);
00267       if (im_on_sub == NULL) {
00268         cpl_msg_error(_id, "sub_image failed\n" );
00269         return -1;
00270       }
00271 
00272 /*
00273     #free memory ------------------------
00274 */
00275       destroy_image(im_on);
00276       destroy_image(im_off);
00277  }
00278  else {
00279       im_on_sub = im_on;
00280  }
00281  
00282 
00283  /*
00284 #---------------------------------------------------------
00285 # convolution with Gaussian if recommended
00286 #---------------------------------------------------------
00287 */
00288 
00289  im_on_gauss = (OneImage*) cpl_calloc (cfg -> nframes, sizeof(OneImage)) ;
00290 
00291 
00292  if (cfg->gaussInd == 1) {
00293    cpl_msg_info(_id,"convolution with Gaussian");
00294       im_on_gauss = convolveNSImageByGauss(im_on_sub, cfg->hw);
00295 
00296       if (im_on_gauss == NULL) {
00297          cpl_msg_error (_id, "could not carry out convolveNSImageByGauss\n" );
00298         return -1;
00299       }
00300       destroy_image(im_on_sub);
00301  }
00302 
00303 
00304      /*
00305 #---------------------------------------------------------
00306 # static bad pixel indication
00307 #---------------------------------------------------------
00308     */
00309  if (cfg->gaussInd == 1) {
00310       im_on_sub = im_on_gauss;
00311  }
00312 
00313  if (cfg->maskInd == 1) {
00314    cpl_msg_info(_id,"static bad pixel indication");
00315       im_mask = load_image(cfg->mask);
00316       if (im_mask == NULL) {
00317         cpl_msg_error (_id, "could not load static bad pixel mask\n" );
00318         return -1;
00319       }
00320       im_on_ind = multImageByMask(im_on_sub, im_mask);
00321       if (im_on_ind == NULL) {
00322         cpl_msg_error (_id, "could not carry out multImageByMask\n" );
00323         return -1;
00324       }
00325       destroy_image(im_mask);
00326       destroy_image(im_on_sub);
00327   
00328  } 
00329  else {
00330       im_on_ind = im_on_sub;
00331  }
00332  cpl_free(im_on_gauss);
00333 
00334 
00335  im_on_ind_imgw=cpl_image_wrap_float(im_on_ind->lx,
00336                                    im_on_ind->ly,
00337                                    im_on_ind->data);
00338  im_on_ind_img=cpl_image_duplicate(im_on_ind_imgw);
00339  cpl_image_unwrap(im_on_ind_imgw);
00340 
00341        if(-1 == sinfoni_pro_dump_ima(im_on_ind_img,raw,sof,cfg->fitsname,
00342            PRO_SLIT,plugin_id,config)) {
00343            cpl_msg_error(_id,"cannot dump ima %s", cfg->fitsname);
00344 
00345          cpl_table_delete(tbl_dist);
00346          cpl_table_delete(qclog_tbl);
00347           cpl_free(list_images);
00348           destroy_image(im_on_ind);
00349           cpl_image_delete(im_on_ind_img);
00350           ns_free (cfg);
00351           cpl_frameset_delete(raw);
00352 
00353            return -1;
00354        }
00355        cpl_image_delete(im_on_ind_img);
00356 
00357 
00358 /*
00359 #---------------------------------------------------------
00360 # do the north - south - test
00361 #---------------------------------------------------------
00362 */
00363        cpl_msg_info(_id,"Do the north - south - test");
00364  /* The following gives seg fault */
00365        /*
00366  distances = cpl_calloc(cfg->nslits, sizeof(float));
00367        */
00368 
00369 
00370 
00371 
00372  distances = north_south_test(im_on_ind, 
00373                              cfg->nslits, 
00374                              cfg->halfWidth, 
00375                              cfg->fwhm , 
00376                              cfg->minDiff, 
00377                              cfg->estimated_dist, 
00378                              cfg->devtol, IMA_PIX_START, IMA_PIX_END );
00379  if (distances == NULL) {
00380    cpl_msg_error(_id,"North South Test distance determination failed");
00381    return -1;
00382  }
00383  /*
00384  parameterToAscii(distances, cfg->nslits - 1, cfg->outName);
00385  */
00386  tbl_dist = cpl_table_new(cfg->nslits - 1);
00387  cpl_table_new_column(tbl_dist,"slitlet_distance", CPL_TYPE_FLOAT);
00388  cpl_table_copy_data_float(tbl_dist,"slitlet_distance",distances);
00389  if(cpl_error_get_code() != CPL_ERROR_NONE) {
00390       cpl_msg_error(_id,"error copying data to table");
00391       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00392       return -1;
00393  }
00394  tbl_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00395  strcpy(tbl_name,cfg->outName);
00396 
00397 
00398 
00399 
00400  qclog_tbl = cpl_table_new(cfg->nslits + 1);
00401  cpl_table_new_column(qclog_tbl,"key_name", CPL_TYPE_STRING);
00402  cpl_table_new_column(qclog_tbl,"key_type", CPL_TYPE_STRING);
00403  cpl_table_new_column(qclog_tbl,"key_value", CPL_TYPE_STRING);
00404  cpl_table_new_column(qclog_tbl,"key_help", CPL_TYPE_STRING);
00405 
00406  qc_dist=cpl_vector_new(cfg->nslits - 1);
00407 
00408  key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00409  key_name  = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00410  for(i=0;i<cfg->nslits - 1;i++) {
00411          sprintf(key_name,"%s%i","QC SL DIST",i);
00412          cpl_table_set_string(qclog_tbl,"key_name",i,key_name);
00413          cpl_table_set_string(qclog_tbl,"key_type",i,"CPL_TYPE_DOUBLE");
00414          sprintf(key_value,"%g",distances[i]);
00415          cpl_table_set_string(qclog_tbl,"key_value",i,key_value);
00416          cpl_table_set_string(qclog_tbl,"key_help",i,"Slitlet distance");
00417 
00418          cpl_vector_set(qc_dist,i,distances[i]);
00419  }
00420  qc_dist_mean=cpl_vector_get_mean(qc_dist);
00421  qc_dist_stdev=cpl_vector_get_stdev(qc_dist);
00422 
00423  cpl_table_set_string(qclog_tbl,"key_name",cfg->nslits-1,"QC SL DISTAVG");
00424  cpl_table_set_string(qclog_tbl,"key_type",cfg->nslits-1,"CPL_TYPE_DOUBLE");
00425  sprintf(key_value,"%g",cpl_vector_get_mean(qc_dist));
00426  cpl_table_set_string(qclog_tbl,"key_value",cfg->nslits-1,key_value);
00427  cpl_table_set_string(qclog_tbl,"key_help",cfg->nslits-1,"Average Slitlet distance");
00428 
00429 
00430 
00431  cpl_table_set_string(qclog_tbl,"key_name",cfg->nslits,"QC SL DISTRMS");
00432  cpl_table_set_string(qclog_tbl,"key_type",cfg->nslits,"CPL_TYPE_DOUBLE");
00433  sprintf(key_value,"%g",qc_dist_stdev);
00434  cpl_table_set_string(qclog_tbl,"key_value",cfg->nslits,key_value);
00435  cpl_table_set_string(qclog_tbl,"key_help",cfg->nslits,"RMS Slitlet distance");
00436 
00437  if(-1 == sinfoni_pro_save_tbl(tbl_dist,raw,sof,tbl_name,
00438            PRO_SLITLETS_DISTANCE,qclog_tbl,plugin_id,config)) {
00439      cpl_msg_error(_id,"cannot dump tbl %s", tbl_name);
00440 
00441          cpl_free(tbl_name);
00442          cpl_vector_delete(qc_dist);
00443          cpl_table_delete(tbl_dist);
00444          cpl_table_delete(qclog_tbl);
00445          cpl_free(key_name);
00446          cpl_free(key_value);
00447      cpl_free(list_images);
00448      cpl_free(distances);
00449      destroy_image(im_on_ind);
00450      ns_free (cfg);
00451          cpl_frameset_delete(raw);
00452           
00453      return -1;
00454  }
00455  cpl_free(tbl_name);
00456  cpl_vector_delete(qc_dist);
00457  cpl_table_delete(tbl_dist);
00458  cpl_table_delete(qclog_tbl);
00459  cpl_free(key_name);
00460  cpl_free(key_value);
00461 /*
00462 # free memory -----------
00463 */
00464  cpl_free(list_images);
00465  cpl_free(distances);
00466  destroy_image(im_on_ind);
00467  ns_free (cfg);
00468  cpl_frameset_delete(raw);
00469  return 0;
00470 }
00471 
00472 
00473 
00474 

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