twiflat.c

00001 /*----------------------------------------------------------------------------
00002    
00003    File name    :       twiflat.c
00004    Author       :   A. Modigliani
00005    Created on   :   Sep 17, 2003
00006    Description  :   Master twiflat creation 
00007 
00008  ---------------------------------------------------------------------------*/
00009 
00010 /*----------------------------------------------------------------------------
00011                                 Includes
00012  ---------------------------------------------------------------------------*/
00013 #include "twiflat.h"
00014 /*----------------------------------------------------------------------------
00015                                 Defines
00016  ---------------------------------------------------------------------------*/
00017 
00018 
00019 /*----------------------------------------------------------------------------
00020                              Function Definitions
00021  ---------------------------------------------------------------------------*/
00022 
00023 /*----------------------------------------------------------------------------
00024    Function     :       twiflat()
00025    In           :       ini_file: file name of according .ini file
00026    Out          :       integer (0 if it worked, -1 if it doesn't) 
00027    Job          :       
00028 
00029 this recipe handles stacks of twilight sky frames, that means it generates 
00030 a spatial and spectral flat field by means of an additional flatfield of an 
00031 integrating sphere. First, a dark frame is subtracted from the twilight frame 
00032 and an off frame from the integrating sphere flat. Then the available images 
00033 are stacked in data cube and the clean mean is taken along the z-axis. The 
00034 static bad pixels are indicated, then a clean mean of each image column is 
00035 taken by leaving the bad pixels. The result is one row of the twilight flat. 
00036 Each column of the integrating sphere flat is normalized by its clean mean to 
00037 get a spectral flat. Finally, each row of this frame is multiplied by the 
00038 twilight flat row.
00039 
00040  ---------------------------------------------------------------------------*/
00041 void change_plist_twiflat (cpl_propertylist*  plist, char* out_name, int n_frames);
00042 
00043 int twiflat (cpl_parameterlist* config,cpl_frameset* sof)
00044 {
00045     const char* _id = "twiflat";
00046     twiflat_config * cfg=NULL ;
00047 
00048     OneImage ** list_twi ;
00049     OneImage ** list_dark ;
00050     OneImage ** list_on ;
00051     OneImage ** list_off ;
00052     OneImage ** im ;
00053 
00054     OneImage * im_twi=NULL ;
00055     OneImage * im_dark=NULL ;
00056     OneImage * im_on=NULL ;
00057     OneImage * im_off=NULL ;
00058     OneImage * im_twi_sub=NULL ;
00059     OneImage * im_twi_bad=NULL ;
00060     OneImage * im_twi_tilt=NULL ;
00061     OneImage * im_on_sub=NULL ;
00062     OneImage * im_on_bad=NULL ;
00063     OneImage * im_on_tilt=NULL ;
00064  
00065     OneImage * im_on_norm ;
00066 
00067     OneImage * flat_field=NULL;
00068     OneImage * mask_im=NULL ;
00069     OneImage* row_twi=NULL;
00070     OneImage* row_on=NULL;
00071 
00072     Vector* row_twi_vec=NULL;
00073     Vector* row_on_vec=NULL;
00074 
00075     OneCube * cube_dark=NULL ;
00076     OneCube * cube_on=NULL ;
00077     OneCube * cube_off=NULL ;
00078     OneCube * cube_twi=NULL ;
00079 
00080 
00081     /* fits_header ** header_list ; */
00082     /* fits_header ** head_list ; */
00083     cpl_propertylist* cplist=NULL;
00084      char* name=NULL;
00085 
00086 
00087     int i = 0;
00088     int typ=0;
00089 
00090     int ntwi=0;
00091     int ndark=0;
00092     int non=0;
00093     int noff=0;
00094  
00095 
00096     cpl_image* ff_img=NULL;
00097     cpl_frameset* raw=NULL;
00098    
00099     /* 
00100        -----------------------------------------------------------------
00101        1) parse the file names and parameters to the twiflat_config data 
00102           structure cfg 
00103        -----------------------------------------------------------------
00104      */
00105 
00106     cfg = parse_cpl_input_twiflat(config,sof,&raw) ;
00107 
00108     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00109       cpl_msg_error(_id,"error parsing cpl input");
00110       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00111       return -1;
00112     }
00113     if (cfg == NULL)
00114     {
00115         cpl_msg_error (_id,"could not parse cpl input!\n") ;
00116         return -1 ;
00117     }
00118 
00119     /* load all the needed parameters */
00120 
00121     if(is_fits_file(cfg->mask)  !=1 ){
00122       cpl_msg_error(_id,"Input file %s cfg->mask is not FITS",cfg->mask);
00123       return -1;
00124     }
00125    
00126 
00127     /*
00128 #---------------------------------------------------------
00129 # Take a clean mean of several images
00130 # input is 1 or more similar images
00131 #---------------------------------------------------------
00132 
00133     */
00134 
00135 
00136     /* allocate memory for lists of object, sky and dithered frames */
00137     list_twi = new_image_list (cfg->ntwi);
00138     if (list_twi == NULL) {
00139       cpl_msg_error(_id,"could not allocate memory\n");
00140       return -1;
00141     }
00142 
00143     list_dark = new_image_list(cfg->ndark);
00144     if (list_dark == NULL) {
00145       cpl_msg_error(_id,"could not allocate memory\n");
00146       return -1;
00147     }
00148     list_on = new_image_list(cfg->non);
00149 
00150     if (list_on == NULL) {
00151       cpl_msg_error(_id,"could not allocate memory\n");
00152       return -1;
00153     }
00154 
00155     list_off = new_image_list(cfg->noff);
00156     if (list_off == NULL) {
00157       cpl_msg_error(_id,"could not allocate memory\n");
00158       return -1;
00159     }
00160 
00161     /* build different image lists for the different cases */
00162     ntwi = 0;
00163     ndark = 0;
00164     non = 0;
00165     noff = 0;
00166     /* head_list = (fits_header**) cpl_calloc(cfg -> nframes, sizeof(fits_header*)); */
00167     /* header_list = (fits_header**) cpl_calloc(cfg -> nframes, sizeof(fits_header*)); */
00168 
00169     im = (OneImage**) cpl_calloc(cfg -> nframes, sizeof(OneImage*));
00170     for (i=0; i<  cfg->nframes; i++) {
00171       name = cfg->framelist[i];
00172       if(is_fits_file(name) != 1) {
00173         cpl_msg_info(_id,"Input file %s is not FITS",name);
00174         return -1;
00175       }
00176       im[i]=load_image( name );
00177     
00178     /* read the fits header to change the gain and noise parameter and others*/
00179       /* head_list[i] = fits_read_header(name); */
00180     }
00181 
00182 
00183     for (i=0; i< cfg->nframes; i++) {
00184       typ = intarray_get_value( cfg->frametype, i );
00185       if (im[i] == NULL) {
00186         cpl_msg_error(_id,"could not load image");
00187         return -1;
00188       }
00189       if (typ == 0) {
00190         insert_image( list_twi, im[i], ntwi );
00191         /* header_list[ntwi] = head_list[i]; */
00192         ntwi = ntwi + 1;
00193       }
00194       else if (typ == 1) {
00195         insert_image( list_dark, im[i], ndark );
00196         ndark = ndark + 1;
00197       }
00198       else if (typ == 2) {
00199         insert_image( list_on, im[i], non );
00200         non = non + 1;
00201       }
00202       else if (typ == 3) {
00203         insert_image( list_off, im[i], noff );
00204         noff = noff + 1;
00205       }
00206       else {
00207         cpl_msg_error(_id,"wrong frame type\n");
00208         return -1;
00209       }
00210     }
00211     if ((cfg->non != non) || 
00212         (cfg->noff != noff) || 
00213         (cfg->ntwi != ntwi) || 
00214         (cfg->ndark != ndark)) {
00215       cpl_msg_error(_id,"something wrong with the number of the different types of frames");
00216       return -1;
00217     }
00218 
00219 
00220     /* create and fill cubes with the different image lists */
00221     cube_twi = list_make_cube(list_twi, ntwi);
00222     if (cube_twi == NULL) {
00223       cpl_msg_error (_id,"could not create data cube!");
00224       return -1;
00225     }
00226 
00227     destroy_list(list_twi);
00228     cube_dark = list_make_cube(list_dark, ndark);
00229 
00230     if (cube_dark == NULL) {
00231       cpl_msg_error (_id,"could not create data cube!");
00232       return -1;
00233     }
00234 
00235     destroy_list(list_dark);
00236     cube_on = list_make_cube(list_on, non);
00237     if (cube_on == NULL) {
00238       cpl_msg_error (_id,"could not create data cube!");
00239       return -1;
00240     }
00241 
00242     destroy_list(list_on);
00243     cube_off = list_make_cube(list_off, noff);
00244     if (cube_off == NULL) {
00245       cpl_msg_error(_id,"could not create data cube!");
00246       return -1;
00247     }
00248 
00249     destroy_list(list_off);
00250     im_twi = average_with_rejection( cube_twi, cfg->loReject, cfg->hiReject );
00251     if (im_twi == NULL) {
00252       cpl_msg_error(_id, "average_with_rejection failed\n" );
00253       return -1;
00254     }
00255 
00256     im_dark = average_with_rejection( cube_dark, 
00257                                       cfg->loReject, 
00258                                       cfg->hiReject );
00259 
00260     if (im_dark == NULL) {
00261       cpl_msg_error (_id, "average_with_rejection failed\n" );
00262       return -1;
00263     }
00264 
00265     im_on = average_with_rejection( cube_on, cfg->loReject, cfg->hiReject );
00266     if (im_on == NULL) {
00267       cpl_msg_error (_id," average_with_rejection failed\n" );
00268       return -1;
00269     }
00270 
00271     im_off = average_with_rejection( cube_off, cfg->loReject, cfg->hiReject );
00272     if (im_off == NULL) {
00273       cpl_msg_error ( _id," average_with_rejection failed\n" );
00274       return -1;
00275     }
00276 
00277     destroy_cube(cube_twi);
00278     destroy_cube(cube_on);
00279     destroy_cube(cube_dark);
00280     destroy_cube(cube_off);
00281 
00282     /* finally, subtract off from on frames and store the result in the 
00283        object cube
00284     */
00285     im_twi_sub = sub_image(im_twi, im_dark);
00286     if (im_twi_sub == NULL) {
00287       cpl_msg_error ( _id," sub_image failed\n" );
00288       return -1;
00289     }
00290 
00291     im_on_sub = sub_image(im_on, im_off);
00292     if (im_on_sub == NULL) {
00293       cpl_msg_error ( _id," sub_image failed\n" );
00294       return -1;
00295     }
00296 
00297     /* free memory */
00298     destroy_image(im_dark);
00299     destroy_image(im_off);
00300     destroy_image(im_twi);
00301     destroy_image(im_on);
00302 
00303     /*
00304 #---------------------------------------------------------
00305 # static bad pixel indication
00306 #---------------------------------------------------------
00307     */
00308 
00309 
00310     mask_im = load_image(cfg->mask);
00311     if (mask_im == NULL) {
00312       cpl_msg_error ( _id," could not load static bad pixel mask\n" );
00313       return -1;
00314     }
00315 
00316     im_twi_bad = multImageByMask(im_twi_sub, mask_im);
00317     if (im_twi_bad == NULL) {
00318       cpl_msg_error ( _id," could not carry out multImageByMask\n" );
00319       return -1;
00320     }
00321 
00322     im_on_bad = multImageByMask(im_on_sub, mask_im);
00323     if (im_on_bad == NULL) {
00324       cpl_msg_error ( _id," could not carry out multImageByMask\n" );
00325       return -1;
00326     }
00327 
00328     destroy_image(mask_im);
00329     destroy_image(im_on_sub);
00330     destroy_image(im_twi_sub);
00331 
00332 
00333     /*
00334 #---------------------------------------------------------
00335 # correction of linear tilt
00336 #---------------------------------------------------------
00337     */
00338 
00339     if (cfg->warpfixInd == 1) {
00340       /*open ASCII file containing the slope parameter and read it */
00341 
00342       im_twi_tilt = image_warp(im_twi_bad, cfg->kernel, cfg->polyFile);
00343       if (im_twi_tilt == NULL) {
00344         cpl_msg_error ( _id," could not carry out image_warp\n" );
00345         return -1;
00346       }
00347 
00348       im_on_tilt = image_warp(im_on_bad, cfg->kernel,cfg->polyFile );
00349       if (im_on_tilt == NULL) {
00350         cpl_msg_error ( _id," could not carry out image_warp\n" );
00351         return -1;
00352       }
00353     
00354       destroy_image(im_twi_bad);
00355       destroy_image(im_on_bad);
00356     }
00357 
00358       /*
00359 #---------------------------------------------------------
00360 # Generation of a spatial and spectral flatfield
00361 #---------------------------------------------------------
00362       */
00363 
00364 
00365     if (cfg->warpfixInd == 1) { 
00366     im_twi = im_twi_tilt;
00367     im_on  = im_on_tilt;
00368     } else {
00369     im_twi = im_twi_bad;
00370     im_on  = im_on_bad;
00371     }
00372 
00373     row_twi = cleanMeanOfColumns(im_twi, cfg->loReject, cfg->hiReject);
00374     if (row_twi == NULL) {
00375     cpl_msg_error ( _id," could not carry out cleanMeanOfColumns\n" );
00376         return -1;
00377     }
00378     row_on = cleanMeanOfColumns(im_on, cfg->loReject, cfg->hiReject);
00379 
00380     if (row_on == NULL) {
00381     cpl_msg_error ( _id," could not carry out cleanMeanOfColumns\n" );
00382     return -1;
00383     }
00384 
00385     row_on_vec = imageToVector(row_on);
00386     row_twi_vec = imageToVector(row_twi);
00387 
00388 
00389       /* normalize the spectra of the halogen lamp flat to remove its 
00390          spatial variations
00391       */
00392 
00393     im_on_norm = divImageByRow( im_on, row_on_vec );
00394     if (im_on_norm == NULL) {
00395     cpl_msg_error ( _id," could not carry out divImageByRow\n" );
00396     return -1;
00397     }
00398 
00399       /* the spatial information of the twilight flat is finally multiplied 
00400           to the halogen lamp spectra
00401       */
00402     flat_field = multRowToImage( im_on_norm, row_twi_vec );
00403 
00404     if (flat_field == NULL) {
00405     cpl_msg_error ( _id," could not carry out divImageByRow\n" );
00406     return -1;
00407     }
00408 
00409 
00410     ff_img=cpl_image_wrap_float(flat_field->lx,flat_field->ly,flat_field->data);
00411 
00412 
00413 
00414 
00415     if(-1 == sinfoni_pro_dump_ima(ff_img,raw,sof,cfg->outName,
00416           PRO_MASTER_TWIFLAT,_id,config)) {
00417       cpl_msg_error(_id,"cannot dump ima %s", cfg->outName);
00418        for (i=0; i<  cfg->nframes; i++) {
00419      /* fits_header_destroy(head_list[i]); */
00420         destroy_image(im[i]);
00421       }
00422 
00423        cplist = cpl_propertylist_new() ;
00424        if ((cpl_error_code)((cplist = cpl_propertylist_load(cfg->outName, 0)) == NULL)) {
00425                cpl_msg_error(_id, "getting header from frame %s",cfg->outName);
00426                cpl_propertylist_delete(cplist) ;
00427                return -1 ;
00428        }
00429 
00430        change_plist_twiflat(cplist, cfg->outName, cfg->nframes); 
00431  
00432        if (cpl_image_save(ff_img, cfg->outName, CPL_BPP_DEFAULT, 
00433                       cplist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00434                cpl_msg_error(_id, "Cannot save the product %s",cfg->outName);
00435                cpl_propertylist_delete(cplist) ;
00436                return -1 ;
00437        }
00438        cpl_propertylist_delete(cplist) ;
00439 
00440 
00441 
00442 
00443        /* free(head_list); */
00444       /* free(header_list); */
00445 
00446       destroy_image(flat_field);
00447       destroy_image(im_on_norm);
00448       destroy_image(im_on);
00449       destroy_image(im_twi);
00450 
00451       destroyVector(row_twi_vec);
00452       destroyVector(row_on_vec);
00453       twiflat_free (cfg);
00454       return -1;
00455     }
00456 
00457     /* free memory */ 
00458     for (i=0; i<  cfg->nframes; i++) {
00459       /* fits_header_destroy(head_list[i]); */
00460         destroy_image(im[i]);
00461     }
00462     /* free(head_list); */
00463     /* free(header_list); */
00464 
00465     destroy_image(flat_field);
00466     destroy_image(im_on_norm);
00467     destroy_image(im_on);
00468     destroy_image(im_twi);
00469 
00470     destroyVector(row_twi_vec);
00471     destroyVector(row_on_vec);
00472     twiflat_free (cfg);
00473 
00474 
00475 
00476     /* destroy_stringarray (cfg->framelist, cfg->nframes); */
00477 
00478     return 0 ;
00479 }
00480 
00481 
00482 void change_plist_twiflat (cpl_propertylist*  plist, char* out_name, int n_frames) {
00483     char* date;
00484     char* text;
00485 
00486     char file_name[FILENAMESZ] ;
00487     char file_tmp[FILENAMESZ] ;
00488     text      = cpl_calloc(MAX_NAME_SIZE,sizeof(char));
00489     date     = get_datetime_iso8601() ;
00490 
00491 
00492     strcpy(file_tmp,out_name);
00493     strcpy(file_name,"twiflat");
00494     strcat(file_name," -f");
00495     /* strcat(file_name,ini_name); */
00496 
00497     if (strlen(file_name) > 30) {
00498        strcpy(file_name,"twiflat");
00499        strcat(file_name," -f ");
00500        /* strcat(file_name,ini_name[strlen(ini_name)-10:strlen(ini_name)]); */ 
00501     }
00502 
00503 
00504 
00505     cpl_propertylist_insert_after_string(plist, "EXPTIME", KEY_NAME_PIPEFILE, out_name) ;
00506     cpl_propertylist_set_comment(plist, KEY_NAME_PIPEFILE,"pipeline result filename" ) ;
00507 
00508     cpl_propertylist_insert_after_string(plist, "EXPTIME", KEY_NAME_HPRO_DID, 
00509             "ESO-VLT-DIC.PRO-1.13") ;
00510     cpl_propertylist_set_comment(plist, KEY_NAME_HPRO_DID,"data dictionary version" ) ;
00511 
00512     cpl_propertylist_insert_after_string(plist, "EXPTIME", KEY_NAME_HPRO_TYPE, 
00513             "REDUCED") ;
00514     cpl_propertylist_set_comment(plist, KEY_NAME_HPRO_TYPE,  "product type" ) ;
00515 
00516     cpl_propertylist_insert_after_string(plist, "EXPTIME", KEY_NAME_HPRO_CATG, 
00517             "MASTER_LAMP_TWIFLAT") ;
00518     cpl_propertylist_set_comment(plist, KEY_NAME_HPRO_CATG,  "product category" ) ;
00519 
00520     cpl_propertylist_insert_after_string(plist, "EXPTIME", KEY_NAME_HPRO_STATUS, 
00521             "OK") ;
00522     cpl_propertylist_set_comment(plist, KEY_NAME_HPRO_STATUS,  "pipeline status" ) ;
00523 
00524     cpl_propertylist_insert_after_string(plist, "EXPTIME", KEY_NAME_HPRO_DATE, 
00525             date) ;
00526     cpl_propertylist_set_comment(plist, KEY_NAME_HPRO_CATG,  "pipeline execution date" ) ;
00527 
00528     cpl_propertylist_insert_after_string(plist, "EXPTIME", KEY_NAME_HPRO_RECID, 
00529             file_name) ;
00530     cpl_propertylist_set_comment(plist, KEY_NAME_HPRO_RECID, "recipe ID" ) ;
00531 
00532     cpl_propertylist_insert_after_string(plist, "EXPTIME", KEY_NAME_HPRO_DRSID, 
00533             PACKAGE_VERSION) ;
00534     cpl_propertylist_set_comment(plist, KEY_NAME_HPRO_DRSID, "data reduction system ID" ) ;
00535    
00536     sprintf(text,"%d",n_frames);
00537     cpl_propertylist_insert_after_string(plist, "EXPTIME", KEY_NAME_HPRO_DATANCOM, 
00538             text) ;
00539     cpl_propertylist_set_comment(plist, KEY_NAME_HPRO_DATANCOM, "# of combined frames" ) ;
00540     
00541     cpl_free(text);
00542 }
00543 

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