sinfo_new_lamp_flats.c

00001 /*
00002  * This file is part of the ESO SINFONI Pipeline
00003  * Copyright (C) 2004,2005 European Southern Observatory
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or
00008  * (at your option) any later version.
00009  *
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00018  */
00019 /*----------------------------------------------------------------------------
00020  
00021      File name    :     sinfo_new_lamp_flats.c
00022    Author       :   A. Modigliani
00023    Created on   :   Sep 29, 2003
00024    Description  : 
00025 
00026  * this step handles stacks of lamp flat fields, 
00027  *  o it takes a clean mean,
00028  *  o subtracts the off- from the on-frames, 
00029  *  o corrects for static bad pixels and normalizes for a master flat field. 
00030  *  o It distinguishes the spectrally dithered frames and 
00031  *  o treats them the same way. 
00032  *  o It can also generate a static bad pixel mask if wished.
00033 
00034 
00035  ---------------------------------------------------------------------------*/
00036 
00037 #ifdef HAVE_CONFIG_H
00038 #  include <config.h>
00039 #endif
00040 
00041 /*----------------------------------------------------------------------------
00042                                 Includes
00043  ---------------------------------------------------------------------------*/
00044 #include "sinfo_new_lamp_flats.h"
00045 #include "sinfo_flat_ini_by_cpl.h"
00046 #include "sinfo_pro_save.h"
00047 #include "sinfo_pro_types.h"
00048 #include "sinfo_functions.h"
00049 #include "sinfo_new_cube_ops.h"
00050 #include "sinfo_error.h"
00051 #include "sinfo_utils_wrappers.h"
00052 #include "sinfo_image_ops.h"
00053 #include "sinfo_utilities.h"
00054 #include "sinfo_globals.h"
00055 /*----------------------------------------------------------------------------
00056                                 Defines
00057  ---------------------------------------------------------------------------*/
00058 
00059 static int 
00060 new_qc_get_cnt(cpl_frameset* on_set, cpl_frameset* of_set, flat_config* cfg);
00061 static int 
00062 new_lamp_flats_det_ncounts(cpl_frameset* raw, flat_config* config);
00063 
00064 static struct {
00065 
00066   double avg_on; 
00067   double std_on;
00068   double avg_of; 
00069   double std_of;
00070   double avg_di; 
00071   double std_di;
00072   double nsat_on;
00073   double noise_on;
00074   double noise_of;
00075   double flux_on;
00076   double    nsat;
00077 
00078 } qc_lampflat;
00079 
00080 /*----------------------------------------------------------------------------
00081                              Function Definitions
00082  ---------------------------------------------------------------------------*/
00090 /*----------------------------------------------------------------------------
00091    @name  sinfo_new_lamp_flats()
00092    @param  plugin_id recipe id
00093    @param  config input parameterlist
00094    @param  sof    input set of frames
00095    @return integer (0 if it worked, -1 if it doesn't) 
00096    @doc 
00097 
00098  * this step handles stacks of lamp flat fields, 
00099  *  o it takes a clean mean,
00100  *  o subtracts the off- from the on-frames, 
00101  *  o corrects for static bad pixels and normalizes for a master flat field. 
00102  *  o It distinguishes the spectrally dithered frames and 
00103  *  o treats them the same way. 
00104  *  o It can also generate a static bad pixel mask if wished.
00105 
00106 
00107  ---------------------------------------------------------------------------*/
00108 
00109 
00110 
00111 
00112 int 
00113 sinfo_new_lamp_flats (const char* plugin_id, 
00114                       cpl_parameterlist* config, cpl_frameset* sof)
00115 {
00116   flat_config * cfg =NULL;
00117   cpl_imagelist* list_object=NULL;
00118   cpl_imagelist* list_dither_object=NULL ;
00119   cpl_imagelist* list_sky=NULL ;
00120   cpl_imagelist* list_dither_sky=NULL;
00121   cpl_image ** im=NULL ;
00122   cpl_image * norm_dith =NULL;
00123   cpl_image * im_obj =NULL;
00124   cpl_image * int_im =NULL;
00125   cpl_image * int_im_dith =NULL;
00126   cpl_image * im_sky =NULL;
00127   cpl_image * im_dither =NULL;
00128   cpl_image * im_obj_sub =NULL;
00129   cpl_image * im_dither_sub =NULL;
00130   cpl_image * im_dither_sky =NULL;
00131   cpl_image ** imMed=NULL ;
00132   cpl_image * colImage =NULL;
00133   cpl_image * mask_im =NULL;
00134   cpl_image * norm =NULL;
00135   cpl_image * compImage =NULL;
00136   cpl_image * maskImage =NULL;
00137   cpl_image * threshIm =NULL;
00138 
00139   char name[MAX_NAME_SIZE];
00140 
00141   int typ;
00142   Stats * stats =NULL;
00143   int i = 0;
00144   int* n=NULL;
00145   int nob =0;
00146   int nsky = 0; 
00147   int nobjdith = 0;
00148   int nskydith = 0; 
00149   int n_badpixels =0;
00150   int pos =0;
00151 
00152   float** slit_edges=NULL;
00153   float local_clean_mean =0.;
00154   float clean_stdev =0.;
00155   float mean_factor =0.;
00156   float val_x=0;
00157   float val_y=0;
00158 
00159   char outNameDither[MAX_NAME_SIZE];
00160   char name_list[MAX_NAME_SIZE];
00161   char tbl_slitpos_name[MAX_NAME_SIZE];
00162 
00163   cpl_table* tbl_slitpos=NULL;
00164   int* status=NULL;
00165 
00166   int n_im_med=0;
00167   cpl_frameset* raw=NULL;
00168 
00169   cpl_table* qclog_tbl=NULL;
00170   int naxis1=0;
00171   int naxis2=0;
00172   double fpn_stdev1=0;
00173   double fpn_stdev2=0;
00174 
00175   int no=0;
00176   float lo_cut=0;
00177   float hi_cut=0;
00178 
00179   /*
00180        -----------------------------------------------------------------
00181        1) parse the file names and parameters to the tilt_config data 
00182           structure cfg 
00183        -----------------------------------------------------------------
00184   */
00185 
00186   cknull_nomsg(raw=cpl_frameset_new());
00187 
00188   cknull(cfg = sinfo_parse_cpl_input_flat(config,sof,&raw),
00189      "could not parse cpl input!");
00190 
00191   if (cfg->interpolInd == 1) {
00192     if(sinfo_is_fits_file(cfg->mask) != 1) {
00193       sinfo_msg_error("Input file %s is not FITS",cfg->mask);
00194       goto cleanup;
00195     }
00196     if (sinfo_is_fits_file(cfg->slitposList) != 1) {
00197       sinfo_msg_error("Input file %s is not FITS",cfg->slitposList);
00198       goto cleanup;
00199     }
00200   }
00201 
00202   /*
00203     #---------------------------------------------------------
00204     # Take a clean mean of several images
00205     # input is 1 or more similar images
00206     #---------------------------------------------------------
00207   */
00208   sinfo_msg("Takes clean mean of several images");
00209   /* #allocate memory for lists of object, sky and dithered frames--*/
00210   cknull(list_object = cpl_imagelist_new (),"could not allocate memory");
00211 
00212   if (cfg->contains_dither == 1) {
00213     cknull(list_dither_object=cpl_imagelist_new(),"could not allocate memory");
00214   }
00215 
00216   if (cfg->contains_sky == 1) {
00217     cknull(list_sky=cpl_imagelist_new(),"could not allocate memory");
00218   }
00219 
00220   if (cfg->contains_dither == 1 && cfg->nditheroff > 0) {
00221     cknull(list_dither_sky=cpl_imagelist_new(),"could not allocate memory");
00222   }
00223 
00224   if (cfg->contains_dither == 0 && cfg->nditheroff > 0){
00225     sinfo_msg_error("please use non-dithered off-frames, remove the 2!");
00226     goto cleanup;
00227   }
00228  
00229 
00230   /* problem with im as image holder: cleanup then does not work */
00231   im = (cpl_image**) cpl_calloc (cfg -> nframes, sizeof(cpl_image*)); 
00232 
00233   for (i=0; i< cfg->nframes; i++) {
00234     strcpy(name,cfg->framelist[i]);
00235     if(sinfo_is_fits_file(name) != 1) {
00236       sinfo_msg_error("PP Input file %s %d is not FITS",name,i);
00237       goto cleanup;
00238     }
00239     im[i]=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00240 
00241   }
00242 
00243   for (i=0; i< cfg->nframes; i++) {
00244     typ = cfg->frametype[i];
00245     pos = cfg->frameposition[i];
00246     cknull(im[i],"could not load image %d",i);
00247     if (pos == 2) {
00248       if (typ == 1) {
00249         cpl_imagelist_set( list_object, cpl_image_duplicate(im[i]), nob );  
00250         nob = nob + 1;
00251       } else {
00252         cpl_imagelist_set( list_sky, cpl_image_duplicate(im[i]), nsky );
00253         nsky = nsky + 1 ;
00254       }
00255     } else {
00256       if (typ == 1) {
00257         cpl_imagelist_set(list_dither_object, 
00258                           cpl_image_duplicate(im[i]), nobjdith );  
00259         nobjdith = nobjdith + 1;
00260       } else {
00261         cpl_imagelist_set( list_dither_sky, 
00262                           cpl_image_duplicate(im[i]), nskydith );
00263         nskydith = nskydith + 1 ;
00264       }
00265     }
00266   }
00267 
00268 
00269   if (nob != cfg->nobj || cfg->noff != nsky || 
00270       nobjdith != cfg->nditherobj || nskydith != cfg->nditheroff) {
00271       sinfo_msg_error("something is wrong with the number of "
00272                       "the different types of frames");
00273       goto cleanup;
00274   }
00275 
00276   /* # create and fill cubes with the different image lists- */
00277   sinfo_msg("Creates and fills cubes with the different image lists");
00278   cknull(list_object,"could not create data cube!");
00279 
00280   if (cfg->contains_dither == 1) {
00281     cknull(list_dither_object,"could not create data cube!");
00282   }
00283   if (cfg->contains_sky == 1 && nsky > 0) {
00284     cknull(list_sky,"could not create data cube!");
00285   }
00286 
00287   if (cfg->contains_dither == 1 && nskydith > 0) {
00288     cknull(list_dither_sky,"could not create data cube!");
00289   }
00290 
00291  
00292   /*-take the average of the different cubes -*/
00293   sinfo_msg("Takes the average of the different cubes");
00294   if (cfg->loReject*cfg->nobj < 1. && cfg->hiReject *cfg->nobj < 1.) {
00295     cknull(im_obj = sinfo_new_average_cube_to_image(list_object ),
00296           "sinfo_averageCubeToImage failed" );
00297   } else {
00298 
00299     no=cpl_imagelist_get_size(list_object);
00300     lo_cut=(floor)(cfg->loReject*no+0.5);
00301     hi_cut=(floor)(cfg->hiReject*no+0.5);
00302     cknull(im_obj=cpl_imagelist_collapse_minmax_create(list_object,
00303                                                       lo_cut,
00304                                                       hi_cut),
00305                  "sinfo_average_with_rejection failed" );
00306   }
00307    sinfo_free_imagelist(&list_object);
00308 
00309   if (cfg->contains_sky == 1) {
00310     if (cfg->loReject * nsky < 1. && cfg->hiReject * nsky < 1.) {
00311       cknull(im_sky = sinfo_new_average_cube_to_image(list_sky ),
00312       "sinfo_new_average_cube_to_image failed" );
00313     } else {
00314 
00315     no=cpl_imagelist_get_size(list_sky);
00316     lo_cut=(floor)(cfg->loReject*no+0.5);
00317     hi_cut=(floor)(cfg->hiReject*no+0.5);
00318     cknull(im_sky=cpl_imagelist_collapse_minmax_create(list_sky,lo_cut,hi_cut),
00319            "sinfo_average_with_rejection failed" );
00320     }
00321     sinfo_free_imagelist(&list_sky);
00322   }
00323 
00324   if (cfg->contains_dither == 1) {
00325     if (cfg->loReject*nobjdith < 1. && cfg->hiReject * nobjdith < 1.) {
00326       cknull(im_dither = sinfo_new_average_cube_to_image(list_dither_object ),
00327              "sinfo_new_average_cube_to_image failed" );
00328    } else {
00329 
00330 
00331      no=cpl_imagelist_get_size(list_dither_object);
00332      lo_cut=(floor)(cfg->loReject*no+0.5);
00333      hi_cut=(floor)(cfg->hiReject*no+0.5);
00334      cknull(im_dither=cpl_imagelist_collapse_minmax_create(list_dither_object,
00335                                                            lo_cut,hi_cut),
00336                                       "sinfo_average_with_rejection failed" );
00337    }
00338    sinfo_free_imagelist(&list_dither_object);
00339   }
00340 
00341   if (cfg->contains_dither == 1 && nskydith > 0 ) {
00342     if (cfg->loReject*nskydith < 1. && cfg->hiReject*nskydith < 1.) {
00343       cknull(im_dither_sky = sinfo_new_average_cube_to_image(list_dither_sky ),
00344             "sinfo_new_average_cube_to_image failed" );
00345     } else {
00346       no=cpl_imagelist_get_size(list_dither_sky);
00347       lo_cut=(floor)(cfg->loReject*no+0.5);
00348       hi_cut=(floor)(cfg->hiReject*no+0.5);
00349       cknull(im_dither_sky=cpl_imagelist_collapse_minmax_create(list_dither_sky,
00350                                                                 lo_cut,hi_cut),
00351                                           "new_average_with_rejection failed" );
00352     }
00353     sinfo_free_imagelist(&list_dither_sky);
00354   }   
00355 
00356   /*
00357   #---------------------------------------------------------
00358   # Subtract the resulting off-frame (sky) from the on-frame 
00359   #-------------------------------------------------------
00360   #finally, subtract off from on frames and store the result in the 
00361   # object cube----------------
00362   */
00363 
00364   sinfo_msg("Subtracts the resulting off-frame (sky) from the on-frame");
00365   if (cfg->contains_sky == 1) {
00366      cknull(im_obj_sub = cpl_image_subtract_create(im_obj, im_sky),
00367             "could not sinfo_sub_image");
00368      sinfo_free_image(&im_obj);
00369      if (((cfg->contains_dither == 1) && (nskydith > 0)) || 
00370          (cfg->contains_dither == 0)) {
00371        sinfo_free_image(&im_sky);
00372      }
00373      im_obj = im_obj_sub;
00374   }
00375 
00376   if (cfg->contains_dither == 1 && nskydith > 0) {
00377      cknull(im_dither_sub=cpl_image_subtract_create(im_dither, im_dither_sky),
00378             "could not sinfo_sub_image");
00379      sinfo_free_image(&im_dither);
00380      sinfo_free_image(&im_dither_sky);
00381      im_dither = im_dither_sub;
00382    } else if (cfg->contains_dither == 1 && 
00383                           nskydith == 0 && 
00384               cfg->contains_sky == 1) {
00385      cknull(im_dither_sub = cpl_image_subtract_create(im_dither, im_sky),
00386         "could not sinfo_sub_image");
00387      sinfo_free_image(&im_dither);
00388      sinfo_free_image(&im_sky);
00389      im_dither = im_dither_sub;
00390    }
00391    /*
00392    #---------------------------------------------------------
00393    # Generating a static bad pixel mask:
00394    # remove the intensity tilt from every column
00395    # and compute the standard deviation on a rectangular zone
00396    #---------------------------------------------------------
00397    */
00398 
00399    sinfo_msg("Generating a static bad pixel mask");
00400    n_im_med = cfg->iterations+1;
00401 
00402    imMed=(cpl_image**) cpl_calloc(n_im_med, sizeof(cpl_image*));
00403 
00404    if (cfg->badInd == 1) {
00405      sinfo_msg("removes the intensity tilt from every column and");
00406      sinfo_msg("computes the standard deviation on a rectangular zone");
00407 
00408      /* this call originates 36 bytes leaks */
00409      cknull(colImage  = sinfo_new_col_tilt( im_obj, cfg->sigmaFactor ),
00410         "sinfo_colTilt failed" );
00411 
00412      cknull(stats = sinfo_new_image_stats_on_rectangle(colImage, 
00413                                    cfg->badLoReject, 
00414                                    cfg->badHiReject, 
00415                                    cfg->llx, 
00416                                    cfg->lly, 
00417                                    cfg->urx, 
00418                                cfg->ury),
00419                                    "sinfo_get_image_stats_on_vig failed\n");
00420      
00421      local_clean_mean = stats->cleanmean;
00422      clean_stdev = stats->cleanstdev;
00423 
00424 
00425      /* indicate pixels with great deviations from the clean mean as bad */
00426      if (cfg->threshInd == 1) {
00427          cknull(threshIm = sinfo_new_thresh_image(colImage, 
00428             local_clean_mean-mean_factor*clean_stdev,
00429             local_clean_mean+mean_factor*clean_stdev),
00430         " sinfo_threshImage failed\n" );
00431      }
00432      if (cfg->threshInd == 0) {
00433         threshIm = colImage;
00434      }
00435 
00436      /*
00437      filter iteratively the images by a sinfo_median filter of the nearest 
00438      neighbors under the condition of a deviation greater than a factor 
00439      times the standard deviation
00440      */
00441 
00442      cknull(imMed[0]= sinfo_new_median_image(threshIm,-cfg->factor*clean_stdev),
00443                                             " sinfo_medianImage failed" );
00444      
00445 
00446        /* AMO check again if here the loop start and ending point are proper */
00447 
00448      for (i=1; i< cfg->iterations+1; i++) {
00449        cknull(imMed[i]=sinfo_new_median_image(imMed[i-1], 
00450                           -cfg->factor*clean_stdev),
00451                                           "sinfo_medianImage failed" );
00452      }
00453 
00454      /* compare the filtered image with the input image */
00455      cknull(compImage=sinfo_new_compare_images(threshIm, 
00456                                                imMed[cfg->iterations], 
00457                                                im_obj),
00458         "sinfo_compareImages failed" );
00459 
00460      /*---generate the bad pixel mask */
00461      n = (int*)cpl_calloc(1,sizeof(int));
00462      cknull(maskImage = sinfo_new_promote_image_to_mask( compImage, n ),
00463             "error in sinfo_promoteImageToMask" );
00464 
00465 
00466      n_badpixels = n[0];
00467      sinfo_msg("No of bad pixels: %d", n_badpixels);
00468 
00469      cknull_nomsg(qclog_tbl = sinfo_qclog_init());
00470      ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC BP-MAP NBADPIX",n_badpixels,
00471                                    "No of bad pixels","%d"));
00472 
00473      ck0(sinfo_pro_save_ima(maskImage,raw,sof,cfg->maskname,
00474                   PRO_BP_MAP,qclog_tbl,plugin_id,config),
00475                "cannot save ima %s", cfg->maskname);
00476 
00477 
00478      /* free memory */
00479      sinfo_free_table(&qclog_tbl);
00480  
00481      sinfo_new_del_Stats(stats);
00482      sinfo_free_int(&n);
00483      sinfo_free_image(&threshIm); /* */
00484      if (cfg->threshInd == 1) {
00485        sinfo_free_image(&colImage);
00486      }
00487      sinfo_free_image(&compImage);
00488      sinfo_free_image(&maskImage);
00489 
00490      for (i=0; i< cfg->iterations+1; i++) {
00491        sinfo_free_image(&imMed[i]);
00492      }
00493 
00494    }
00495 
00496    cpl_free(imMed);
00497    imMed=NULL;
00498 
00499    /*
00500    #---------------------------------------------------------
00501    # Master flat field: static bad pixel correction and normalizing
00502    #---------------------------------------------------------
00503    */
00504 
00505    sinfo_msg("Creates a Master flat field");
00506    if (cfg->interpolInd == 1) {
00507      cknull(mask_im = cpl_image_load(cfg->mask,CPL_TYPE_FLOAT,0,0),
00508         "could not load static bad pixel mask" );
00509 
00510       /* open the ASCII list of the slitlet positions */
00511 
00512       /* slit_edges = sinfo_new_2Dfloat_array(32, 2) ; */
00513       slit_edges = (float **) cpl_calloc( 32, sizeof (float*) ) ;
00514       for ( i = 0 ; i < 32 ; i++ )
00515         {
00516             slit_edges[i] = (float *) cpl_calloc( 2, sizeof (float)) ;
00517         }
00518       /*READ TFITS TABLE*/
00519       if(sinfo_is_fits_file(cfg->slitposList) !=1 ) {
00520           sinfo_msg_error("Input file %s is not FITS", cfg->slitposList);
00521           goto cleanup;
00522       }
00523       strcpy(tbl_slitpos_name,cfg->slitposList);
00524       check(tbl_slitpos = cpl_table_load(tbl_slitpos_name,1,0),
00525          "error loading tbl %s",tbl_slitpos_name);
00526 
00527       for (i =0 ; i< 32; i++){
00528             val_x=cpl_table_get_double(tbl_slitpos,"pos1",i,status);
00529             val_y=cpl_table_get_double(tbl_slitpos,"pos2",i,status);
00530             slit_edges[i][0]=val_x;
00531         slit_edges[i][1]=val_y;
00532       }
00533       sinfo_free_table(&tbl_slitpos);
00534 
00535       cknull(int_im = sinfo_interpol_source_image (im_obj, mask_im, 
00536                                                  cfg->maxRad, slit_edges),
00537                          "could not carry out sinfo_interpolSourceImage" );
00538       
00539       sinfo_free_image(&im_obj);
00540       cknull(norm = sinfo_new_normalize_to_central_pixel(int_im),
00541              "could not normalize flatfield" );
00542       sinfo_free_image(&int_im);
00543       im_obj = norm; 
00544 
00545       if (cfg->contains_dither == 1) {
00546         cknull(int_im_dith = sinfo_interpol_source_image(im_dither, 
00547                                mask_im,
00548                                cfg->maxRad,
00549                                slit_edges),
00550            "could not carry out sinfo_interpolSourceImage" );
00551       
00552         cpl_image_delete(im_dither);
00553         cknull(norm_dith = sinfo_new_normalize_to_central_pixel(int_im_dith),
00554            "could not normalize flatfield" );
00555         sinfo_free_image(&int_im_dith);
00556         im_dither = norm_dith;
00557       }
00558       /* sinfo_new_destroy_2Dfloatarray(slit_edges, 32); */    
00559       for ( i = 0 ; i < 32 ; i++ )
00560         {
00561             cpl_free( slit_edges[i] );
00562         }
00563       cpl_free( slit_edges ) ;
00564       sinfo_free_image(&mask_im);
00565 
00566    }
00567 
00568    if (cfg->interpolInd != 1) {
00569      cknull(norm = sinfo_new_normalize_to_central_pixel(im_obj),
00570             "could not normalize flatfield" );
00571       sinfo_free_image(&im_obj);
00572       im_obj = norm; 
00573 
00574       if (cfg->contains_dither == 1) {
00575         cknull(norm_dith = sinfo_new_normalize_to_central_pixel(im_dither),
00576            "could not normalize flatfield" );
00577         sinfo_free_image(&im_dither);
00578         im_dither = norm_dith;
00579       }
00580    }
00581 
00582    naxis1=cpl_image_get_size_x(im_obj);
00583    naxis2=cpl_image_get_size_y(im_obj);
00584 
00585 
00586    if(cfg->qc_fpn_xmin1 < 1) {
00587      sinfo_msg_error("qc_ron_xmin < 1");
00588      goto cleanup;
00589    }
00590 
00591    if(cfg->qc_fpn_xmax1 > naxis1) {
00592      sinfo_msg_error("qc_ron_xmax < %d",naxis1);
00593      goto cleanup;
00594    }
00595 
00596    if(cfg->qc_fpn_ymin1 < 1) {
00597      sinfo_msg_error("qc_ron_ymin < 1");
00598      goto cleanup;
00599    }
00600 
00601    if(cfg->qc_fpn_ymax1 > naxis2) {
00602      sinfo_msg_error("qc_ron_ymax < %d",naxis2);
00603      goto cleanup;
00604    }
00605    fpn_stdev1 = cpl_image_get_stdev_window(im_obj,
00606                        cfg->qc_fpn_xmin1,
00607                        cfg->qc_fpn_ymin1,
00608                        cfg->qc_fpn_xmax1,
00609                        cfg->qc_fpn_ymax1);
00610 
00611 
00612    if(cfg->qc_fpn_xmin2 < 1) {
00613      sinfo_msg_error("qc_ron_xmin < %d",1);
00614      goto cleanup;
00615    }
00616 
00617 
00618    if(cfg->qc_fpn_xmax2 > naxis1) {
00619      sinfo_msg_error("qc_ron_xmax < %d",naxis1);
00620      goto cleanup;
00621    }
00622 
00623    if(cfg->qc_fpn_ymin2 < 1) {
00624      sinfo_msg_error("qc_ron_ymin < 1");
00625      goto cleanup;
00626    }
00627 
00628    if(cfg->qc_fpn_ymax2 > naxis2) {
00629      sinfo_msg_error("qc_ron_ymax < %d",naxis2);
00630      goto cleanup;
00631    }
00632    fpn_stdev2 = cpl_image_get_stdev_window(im_obj,
00633                        cfg->qc_fpn_xmin2,
00634                        cfg->qc_fpn_ymin2,
00635                        cfg->qc_fpn_xmax2,
00636                        cfg->qc_fpn_ymax2);
00637 
00638 
00639 
00640    ck0(new_lamp_flats_det_ncounts(raw,cfg),"error computing number of counts");
00641    cknull_nomsg(qclog_tbl = sinfo_qclog_init());
00642    ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC SPECFLAT NCNTSAVG",
00643                                     qc_lampflat.avg_di,"Average counts","%g"));
00644    ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC SPECFLAT NCNTSSTD",
00645                                     qc_lampflat.std_di,"Stdev counts","%g"));
00646    ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC SPECFLAT OFFFLUX",
00647                                     qc_lampflat.avg_of,
00648                                     "Average flux off frames","%g"));
00649 
00650     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,
00651                                      "QC LFLAT FPN1",
00652                                      fpn_stdev1,
00653                                      "Fixed Pattern Noise of combined frames",
00654                                      "%f"));
00655 
00656     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,
00657                                      "QC LFLAT FPN2",
00658                                      fpn_stdev2,
00659                                      "Fixed Pattern Noise of combined frames",
00660                                      "%f"));
00661 
00662     ck0(sinfo_pro_save_ima(im_obj,raw,sof,cfg->outName,
00663                PRO_MASTER_FLAT_LAMP,qclog_tbl,plugin_id,config),
00664                            "cannot save ima %s", cfg->outName);
00665 
00666 
00667     sinfo_free_table(&qclog_tbl);
00668     sinfo_free_image(&im_obj);
00669 
00670  
00671     if (cfg->contains_dither == 1) {
00672 
00673       if (strstr(cfg->outName, ".fits" ) != NULL ) {
00674 
00675     snprintf(name_list, MAX_NAME_SIZE-1,"%s%s", sinfo_new_get_rootname(cfg->outName),
00676                 "_dither");
00677         strcpy(outNameDither,name_list);
00678         strcat(outNameDither,strstr(cfg->outName,".fits"));
00679        
00680       } else {
00681         strcpy(outNameDither,cfg->outName);
00682         strcat(outNameDither,"_dither");
00683       }
00684 
00685 
00686       naxis1=cpl_image_get_size_x(im_dither);
00687       naxis2=cpl_image_get_size_y(im_dither);
00688 
00689 
00690       if(cfg->qc_fpn_xmin1 < 1) {
00691      sinfo_msg_error("qc_ron_xmin1 < 1");
00692          goto cleanup;
00693       }
00694 
00695       if(cfg->qc_fpn_xmax1 > naxis1) {
00696      sinfo_msg_error("qc_ron_xmax1 < %d",naxis1);
00697          goto cleanup;
00698       }
00699 
00700       if(cfg->qc_fpn_ymin1 < 1) {
00701      sinfo_msg_error("qc_ron_ymin1 < 1");
00702          goto cleanup;
00703       }
00704 
00705       if(cfg->qc_fpn_ymax1 > naxis2) {
00706      sinfo_msg_error("qc_ron_ymax1 < %d",naxis2);
00707          goto cleanup;
00708       }
00709 
00710  
00711       fpn_stdev1 = cpl_image_get_stdev_window(im_dither,
00712                           cfg->qc_fpn_xmin1,
00713                           cfg->qc_fpn_ymin1,
00714                           cfg->qc_fpn_xmax1,
00715                           cfg->qc_fpn_ymax1);
00716 
00717       if(cfg->qc_fpn_xmin2 < 1) {
00718      sinfo_msg_error("qc_ron_xmin2 < 1");
00719          goto cleanup;
00720       }
00721 
00722       if(cfg->qc_fpn_xmax2 > naxis1) {
00723      sinfo_msg_error("qc_ron_xmax2 < %d",naxis1);
00724          goto cleanup;
00725       }
00726 
00727       if(cfg->qc_fpn_ymin2 < 1) {
00728      sinfo_msg_error("qc_ron_ymin2 < 1");
00729          goto cleanup;
00730       }
00731 
00732       if(cfg->qc_fpn_ymax2 > naxis2) {
00733      sinfo_msg_error("qc_ron_ymax2 < %d",naxis2);
00734          goto cleanup;
00735       }
00736  
00737       fpn_stdev2 = cpl_image_get_stdev_window(im_dither,
00738                           cfg->qc_fpn_xmin2,
00739                           cfg->qc_fpn_ymin2,
00740                           cfg->qc_fpn_xmax2,
00741                           cfg->qc_fpn_ymax2);
00742 
00743 
00744       ck0(new_lamp_flats_det_ncounts(raw,cfg),"error computing ncounts");
00745       cknull_nomsg(qclog_tbl = sinfo_qclog_init());
00746       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC SPECFLAT NCNTSAVG",
00747                                     qc_lampflat.avg_di,"Average counts","%g"));
00748 
00749       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC SPECFLAT NCNTSSTD",
00750                                     qc_lampflat.std_di,"Stdev counts","%g"));
00751 
00752       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC SPECFLAT OFFFLUX",
00753                            qc_lampflat.avg_of,"Average flux off frames","%g"));
00754 
00755       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC LFLAT FPN1",fpn_stdev1,
00756                                "Fixed Pattern Noise of combined frames","%f"));
00757 
00758       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC LFLAT FPN2",fpn_stdev2,
00759                                "Fixed Pattern Noise of combined frames","%f"));
00760  
00761 
00762       ck0(sinfo_pro_save_ima(im_dither,raw,sof,outNameDither,
00763                  PRO_MASTER_FLAT_LAMP,qclog_tbl,plugin_id,config),
00764                     "cannot save ima %s", outNameDither);
00765 
00766 
00767       sinfo_free_table(&qclog_tbl);
00768       sinfo_free_image(&im_dither);
00769 
00770     }
00771 
00772 
00773     /* could be done earlier? */
00774     sinfo_free_image_array(&im,cfg->nframes);
00775     sinfo_free_frameset(&raw);
00776     sinfo_flat_free(&cfg);
00777     return 0;
00778 
00779 cleanup:
00780 
00781 
00782     /* free memory */
00783     if(slit_edges != NULL) {
00784       for ( i = 0 ; i < 32 ; i++ )
00785         {
00786       if(slit_edges[i] != NULL) {
00787             cpl_free( slit_edges[i] );
00788       }
00789           slit_edges[i]=NULL;
00790         }
00791       cpl_free( slit_edges ) ;
00792     }
00793     sinfo_free_image(&mask_im);
00794     sinfo_free_table(&qclog_tbl);
00795     if(stats) {
00796       sinfo_new_del_Stats(stats);
00797     }
00798     sinfo_free_int(&n);
00799     sinfo_free_image(&threshIm);
00800     sinfo_free_image(&maskImage);
00801     if(imMed != NULL) sinfo_free_image_array(&imMed,cfg->iterations);
00802     if(stats!= NULL) {
00803         sinfo_new_del_Stats(stats);
00804         stats=NULL;
00805     }
00806     sinfo_free_image(&compImage);
00807     sinfo_free_image(&colImage);
00808     sinfo_free_imagelist(&list_dither_object);
00809     sinfo_free_imagelist(&list_dither_sky);
00810     if(list_sky != NULL) {
00811        sinfo_free_imagelist(&list_sky);
00812     }
00813     sinfo_free_imagelist(&list_object);
00814     sinfo_free_image(&im_dither);
00815     sinfo_free_image(&im_dither_sky);
00816     sinfo_free_image(&im_obj);
00817     sinfo_free_image(&im_sky);
00818     if(im != NULL) sinfo_free_image_array(&im,cfg->nframes);
00819     sinfo_free_frameset(&raw);
00820     if(cfg != NULL) {
00821       sinfo_flat_free(&cfg);
00822     }
00823     return -1;
00824 
00825 }
00826 
00827 static int
00828 new_lamp_flats_det_ncounts(cpl_frameset* raw, flat_config* cfg)
00829 {
00830  int i=0;
00831  int j=0;
00832 
00833  int nraw=0;
00834  int non=0;
00835  int noff=0;
00836 
00837  double mjd_on=0;
00838  double mjd_of=0;
00839  double mjd_of_frm=0;
00840 
00841  char filename[MAX_NAME_SIZE];
00842 
00843  cpl_frame* frm=NULL;
00844  cpl_frame* frm_dup=NULL;
00845  cpl_frame* on_frm=NULL;
00846  cpl_frame* of_frm=NULL;
00847  cpl_frame* tmp_of_frm=NULL;
00848 
00849 
00850  cpl_frameset* on_set=NULL;
00851  cpl_frameset* of_set=NULL;
00852  cpl_frameset* wrk_set=NULL;
00853 
00854  on_set=cpl_frameset_new();
00855  of_set=cpl_frameset_new();
00856 
00857  nraw = cpl_frameset_get_size(raw);
00858 
00859  for (i=0; i< nraw; i++) {
00860    frm = cpl_frameset_get_frame(raw,i);
00861    frm_dup = cpl_frame_duplicate(frm);
00862    if(sinfo_frame_is_on(frm) == 1) {
00863      cpl_frameset_insert(on_set,frm_dup);
00864      non++;
00865    } else {
00866      cpl_frameset_insert(of_set,frm_dup);
00867      noff++;
00868    }
00869  }
00870 
00871 
00872  if (non == noff) {
00873    new_qc_get_cnt(on_set,of_set,cfg);
00874 
00875  } else if (non == 0) {
00876    sinfo_msg("non == 0");
00877    sinfo_msg_warning("QC SPECFLAT NCNTAVG/NCTNTSTD/OFFFLUX=0 ");
00878  
00879  } else if ( noff == 0 ) {
00880    sinfo_msg("noff == 0");
00881    sinfo_msg_warning("QC SPECFLAT NCNTAVG/NCTNTSTD/OFFFLUX=0 ");
00882 
00883  } else {
00884 
00885    sinfo_msg_warning("non != noff, => QC SPECFLAT NCNTAVG/NCTNTSTD/OFFFLUX=0 ");
00886 
00887     for (i=0;i<non;i++) {
00888        wrk_set=cpl_frameset_new();
00889        on_frm=cpl_frameset_get_frame(on_set,i);
00890        mjd_on=sinfo_get_mjd_obs(on_frm);
00891        of_frm=cpl_frameset_get_frame(of_set,0);
00892        mjd_of=sinfo_get_mjd_obs(of_frm);
00893        strcpy(filename,cpl_frame_get_filename(of_frm));
00894        for (j=1;j<noff;j++) {
00895           tmp_of_frm = cpl_frameset_get_frame(of_set,j);
00896           mjd_of_frm = sinfo_get_mjd_obs(tmp_of_frm);
00897 
00898           if(1000.*(mjd_of_frm-mjd_on)*(mjd_of_frm-mjd_on) <
00899              1000.*(mjd_of-    mjd_on)*(mjd_of-    mjd_on) ) {
00900          mjd_of=mjd_of_frm;
00901              of_frm=cpl_frame_duplicate(tmp_of_frm);
00902       }
00903        }
00904        strcpy(filename,cpl_frame_get_filename(of_frm));
00905        frm_dup=cpl_frame_duplicate(of_frm);
00906        cpl_frameset_insert(wrk_set,frm_dup);
00907        strcpy(filename,cpl_frame_get_filename(of_frm));
00908     }
00909     /* Commented out as algorithm non robust if non != noff */
00910     new_qc_get_cnt(on_set,wrk_set,cfg); 
00911 
00912  }
00913 
00914  cpl_frameset_delete(wrk_set);
00915  cpl_frameset_delete(on_set);
00916  cpl_frameset_delete(of_set);
00917  return 0;
00918 
00919 
00920 }
00921 
00922 static int
00923 new_qc_get_cnt(cpl_frameset* on_set, cpl_frameset* of_set, flat_config* cfg)
00924 {
00925 
00926   int i=0;
00927   int nsat=0;
00928   int non=0;
00929   int nof=0;
00930   int nfr=0;
00931 
00932   char name[MAX_NAME_SIZE];
00933   cpl_vector* vec_on=NULL;
00934   cpl_vector* vec_of=NULL;
00935   cpl_vector* vec_di=NULL;
00936   cpl_vector* vec_nsat=NULL;
00937   cpl_frame* on_frm=NULL;
00938   cpl_frame* of_frm=NULL;
00939 
00940   cpl_image* dif_ima=NULL;
00941   cpl_image* on_ima=NULL;
00942   cpl_image* of_ima=NULL;
00943   cpl_image* tmp_ima=NULL;
00944 
00945   double med=0;
00946     non = cpl_frameset_get_size(on_set);
00947     nof = cpl_frameset_get_size(of_set);
00948     nfr = (non <= nof) ? non : nof;
00949     vec_on = cpl_vector_new(nfr);
00950     vec_of = cpl_vector_new(nfr);
00951     vec_di = cpl_vector_new(nfr);
00952     vec_nsat = cpl_vector_new(nfr);
00953 
00954 
00955     for (i=0; i< nfr; i++) {
00956       on_frm = cpl_frameset_get_frame(on_set,i);
00957       strcpy(name,cpl_frame_get_filename(on_frm));
00958       on_ima = cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00959       med= cpl_image_get_median(on_ima);
00960       cpl_vector_set(vec_on,i,med);
00961 
00962       tmp_ima = cpl_image_duplicate(on_ima);
00963       cpl_image_threshold(tmp_ima,CPL_PIXEL_MINVAL,
00964                                cfg->qc_thresh_max,0,1);
00965       nsat=cpl_image_get_flux(tmp_ima);
00966       cpl_vector_set(vec_nsat,i,nsat);
00967 
00968    /* Are you sure to have same frames off as on ? */
00969       of_frm = cpl_frameset_get_frame(of_set,i);
00970       strcpy(name,cpl_frame_get_filename(of_frm));
00971       of_ima = cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00972       med= cpl_image_get_median(of_ima);
00973       cpl_vector_set(vec_of,i,med);
00974       dif_ima = cpl_image_subtract_create(on_ima,of_ima);
00975       med= cpl_image_get_median(dif_ima);
00976       cpl_vector_set(vec_di,i,med);
00977 
00978     cpl_image_delete(on_ima);
00979     cpl_image_delete(of_ima);
00980     cpl_image_delete(dif_ima);
00981     cpl_image_delete(tmp_ima);
00982     }
00983     qc_lampflat.avg_on=cpl_vector_get_mean(vec_on);
00984     qc_lampflat.avg_of=cpl_vector_get_mean(vec_of);
00985     qc_lampflat.avg_di=cpl_vector_get_mean(vec_di);
00986     if(nfr > 1 ) {
00987        qc_lampflat.std_on=cpl_vector_get_stdev(vec_on);
00988        qc_lampflat.std_of=cpl_vector_get_stdev(vec_of);
00989        qc_lampflat.std_di=cpl_vector_get_stdev(vec_di);
00990     }
00991     qc_lampflat.nsat=cpl_vector_get_mean(vec_nsat);
00992     cpl_vector_delete(vec_on);
00993     cpl_vector_delete(vec_of);
00994     cpl_vector_delete(vec_di);
00995     cpl_vector_delete(vec_nsat);
00996     /*
00997     sinfo_msg( "sinfo_qc_get_cnt","avg_on=%g std_on=%g ",
00998                       qc_lampflat.avg_on,qc_lampflat.std_on);
00999     sinfo_msg( "sinfo_qc_get_cnt","avg_of=%g std_of=%g ",
01000                       qc_lampflat.avg_of,qc_lampflat.std_of);
01001     sinfo_msg( "sinfo_qc_get_cnt","avg_di=%g std_di=%g ",
01002                       qc_lampflat.avg_di,qc_lampflat.std_di);
01003     sinfo_msg( "sinfo_qc_get_cnt","nsat=%g ",qc_lampflat.nsat);
01004     */
01005     return 0;
01006 }
01007 

Generated on Wed Jan 17 08:33:43 2007 for SINFONI Pipeline Reference Manual by  doxygen 1.4.4