sinfo_new_cubes_coadd.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_cubes_coadd.c
00022    Author       :   J. Schreiber
00023    Created on   :   December 3, 2003
00024    Description  :   Creates data cubes or merges data cubes 
00025                         out of jittered object-sky
00026                         nodding observations 
00027  ---------------------------------------------------------------------------*/
00028 #ifdef HAVE_CONFIG_H
00029 #  include <config.h>
00030 #endif
00031 
00032 /*----------------------------------------------------------------------------
00033                                 Includes
00034  ---------------------------------------------------------------------------*/
00035 #include "sinfo_new_cubes_coadd.h"
00036 #include "sinfo_pfits.h"
00037 #include "sinfo_pro_save.h"
00038 #include "sinfo_objnod_ini_by_cpl.h"
00039 #include "sinfo_functions.h"
00040 #include "sinfo_utilities_scired.h" 
00041 #include "sinfo_wave_calibration.h"
00042 #include "sinfo_cube_construct.h"
00043 #include "sinfo_error.h"
00044 #include "sinfo_utils_wrappers.h"
00045 
00046 /*----------------------------------------------------------------------------
00047                                 Defines
00048  ---------------------------------------------------------------------------*/
00049 #define PI_NUMB        (3.1415926535897932384626433832795) /* pi */
00050 
00051 
00052 /*----------------------------------------------------------------------------
00053                              Function Definitions
00054  ---------------------------------------------------------------------------*/
00055 
00065 /*----------------------------------------------------------------------------
00066    Function     : sinfo_new_cubes_coadd()
00067    In           : ini_file: file name of according .ini file
00068    Out          : integer (0 if it worked, -1 if it doesn't) 
00069    Job          : this routine carries through the data cube creation of an 
00070                   object science observation using object-sky nodding
00071                   and jittering. This script expects jittered frames that
00072           were already sky-subtracted
00073                   averaged, flatfielded, spectral tilt corrected and 
00074         interleaved if necessary
00075  ---------------------------------------------------------------------------*/
00076 int 
00077 sinfo_new_cubes_coadd (const char* plugin_id,cpl_parameterlist* config, 
00078             cpl_frameset* sof, const char* procatg)
00079 {
00080 
00081     object_config * cfg=NULL ;
00082     cpl_image * im=NULL ;
00083     cpl_image* jitter_image=NULL;
00084     cpl_imagelist  ** cube_tmp=NULL ;
00085 
00086     cpl_imagelist  ** cubeobject=NULL ;
00087     cpl_imagelist  * jittercube=NULL ;
00088     cpl_imagelist  * maskcube=NULL ;
00089     cpl_propertylist* plist=NULL;
00090     int sky_cor=0;
00091     int pdensity=0;
00092 
00093     int i=0;
00094     int n=0 ;
00095     int partind = 0 ;
00096     int centralpix=0 ;
00097     int z_siz=0;
00098     int z_min=0;
00099     int z_max=0;
00100     int z=0;
00101     int z_stp=100;
00102     int scales_sky=0;
00103     int ks_clip=0;
00104     double kappa=2.0;
00105 
00106     float ref_offx=0;
00107     float ref_offy=0;
00108     float center_x=0;
00109     float newcenter_x=0 ;
00110     float center_y=0;
00111     float newcenter_y=0 ;
00112     float cd1_1=0;
00113     float cd1_2=0;
00114     float cd2_1=0;
00115     float cd2_2=0;
00116     double pixelscale=0;
00117     double angle=0;
00118     float radangle=0 ;
00119     double exptime=0 ;
00120 
00121     double * times=NULL ;
00122     float * offsetx=NULL;
00123     float * offsety=NULL;
00124 
00125     double dis=0 ;
00126     double centralLambda=0 ;
00127 
00128     char name_jitter[MAX_NAME_SIZE] ;
00129     char pro_mjit[MAX_NAME_SIZE];
00130     char pro_obs[MAX_NAME_SIZE];
00131     char pro_med[MAX_NAME_SIZE];
00132 
00133 
00134     char * name=NULL ;
00135     char * partname=NULL;
00136     char * partname2=NULL ;
00137     char file_name[MAX_NAME_SIZE];
00138     int vllx=0;
00139     int vlly=0;
00140     int vurx=0;
00141     int vury=0;
00142 
00143     int onp=0;
00144     int j=0;
00145     cpl_image* j_img=NULL;
00146     cpl_image* m_img=NULL;
00147     cpl_table* qclog_tbl=NULL;
00148     cpl_image* ill_cor=NULL;
00149     cpl_frame* frame=NULL;
00150     cpl_frameset* stk=NULL;
00151     cpl_parameter* p=NULL;
00152 
00153      if (strcmp(procatg,PRO_COADD_STD) == 0) {
00154     strcpy(pro_mjit,PRO_MASK_COADD_STD);
00155     strcpy(pro_obs,PRO_OBS_STD);
00156     strcpy(pro_med,PRO_MED_COADD_STD);
00157 
00158     } else if (strcmp(procatg,PRO_COADD_PSF) == 0) {
00159     strcpy(pro_mjit,PRO_MASK_COADD_PSF);
00160     strcpy(pro_obs,PRO_OBS_PSF);
00161     strcpy(pro_med,PRO_MED_COADD_PSF);
00162     } else {
00163     strcpy(pro_mjit,PRO_MASK_COADD_OBJ);
00164     strcpy(pro_obs,PRO_OBS_OBJ);
00165     strcpy(pro_med,PRO_MED_COADD_OBJ);
00166     }
00167 
00168     /*----parse input data and parameters to set cube_config cfg---*/
00169     check_nomsg(stk = cpl_frameset_new());
00170 
00171     cknull(cfg = sinfo_parse_cpl_input_objnod(config,sof,&stk),
00172        "Error setting parameter configuration");
00173  
00174     ck0(sinfo_check_input_data(cfg),"error checking input");
00175 
00176     if ( cfg->jitterind == 1 )
00177     {
00178         cknull(times = (double*) cpl_calloc (cfg->nframes, sizeof (double)), 
00179            " could not allocate memory!") ;
00180  
00181         cknull(offsetx = (float*) cpl_calloc (cfg->nframes, sizeof(float)),
00182            " could not allocate memory!") ;
00183 
00184         cknull(offsety = (float*) cpl_calloc (cfg->nframes, sizeof(float)),
00185            " could not allocate memory!") ;
00186     }
00187   
00188     if (cfg->jitterind == 0)
00189     {
00190         if ( NULL != (partname = strtok(cfg->outName, "."))) 
00191         {
00192             partname2 = strtok (NULL, ".") ;
00193             partind = 1 ;
00194         }
00195     }
00196 
00197 
00198     ck0(sinfo_auto_size_cube2(cfg,&ref_offx,&ref_offy),"Error resizing cube");
00199  
00200     if(NULL != cpl_frameset_find(sof,PRO_ILL_COR)) {
00201      frame = cpl_frameset_find(sof,PRO_ILL_COR);
00202      ill_cor=cpl_image_load(cpl_frame_get_filename(frame),CPL_TYPE_FLOAT,0,0);
00203     } else {
00204      sinfo_msg("Illumination correction image not provided");
00205      cpl_error_reset();
00206     }
00207  
00208      for ( n = 0 ; n < cfg->nframes ; n++ )
00209      {
00210 
00211         sinfo_msg_debug("Read FITS information");
00212         name = cfg->framelist[n] ;    
00213         if (n == 0)
00214     {
00215         strcpy (name_jitter, name) ;
00216     }   
00217         if( sinfo_is_fits_file(name) != 1) {
00218           sinfo_msg_error("Input file %s is not FITS",name);
00219           goto cleanup;
00220     }
00221 
00222     
00223     /* get some header values and compute the CD-sinfo_matrix */
00224         plist=cpl_propertylist_load(name,0);
00225     pixelscale = sinfo_pfits_get_pixscale(plist) /2;
00226     angle = sinfo_pfits_get_posangle(plist) ;
00227         /* in PUPIL data there is not posangle info: we reset the error */
00228         if(cpl_error_get_code() != CPL_ERROR_NONE) {
00229           cpl_error_reset();
00230         }
00231 
00232         sinfo_free_propertylist(&plist);
00233     radangle = angle * PI_NUMB / 180. ;
00234     cd1_1 = cos(radangle) ;
00235     cd1_2 = sin(radangle) ;
00236     cd2_1 = -sin(radangle) ;
00237     cd2_2 = cos(radangle) ;
00238     
00239     sinfo_msg_debug("frame no.: %d, name: %s\n", n, name) ;
00240     cknull(im = cpl_image_load(name,CPL_TYPE_FLOAT,0,0),
00241            " could not load frame %s!",name) ;
00242 
00243     if (cfg->jitterind == 1)
00244     {
00245       exptime = sinfo_pfits_get_ditndit(name) ;
00246 
00247       if (exptime == FLAG)
00248         {
00249           sinfo_msg_error("could not read fits header keyword exptime!");
00250           goto cleanup;
00251         }
00252       times[n] = exptime ;
00253 
00254           ck0(sinfo_new_assign_offset(n,name,offsetx,offsety,
00255                                       ref_offx,ref_offy),
00256                                       "Error assigning offsets");
00257 
00258     }
00259         sinfo_free_image(&im);
00260 
00261     } /* end loop over n (nframes) */
00262  
00263 
00264     /* leak free */
00265     if(cfg->jitterind == 0) {
00266       goto exit;
00267     }   
00268      
00269     /* Here in case of autojitter we estimate the sky */
00270     if( (cfg->size_x*cfg->size_y*cfg->nframes) > (130*130*21) ) {
00271       sinfo_msg_warning("Coadd cube size:%d,%d. N frames: %d",
00272                            cfg->size_x,cfg->size_y,cfg->nframes);
00273       sinfo_msg_warning("Max allowed should be such "
00274                         "that sixeX*sixeY*Nframes < 130*130*21");
00275       goto exit;
00276     } 
00277 
00278     if ( cfg->jitterind == 1 )
00279     {  
00280         check_nomsg(p = cpl_parameterlist_find(config, "sinfoni.objnod.vllx"));
00281         check_nomsg(vllx = cpl_parameter_get_int(p));
00282         check_nomsg(p = cpl_parameterlist_find(config, "sinfoni.objnod.vlly"));
00283         check_nomsg(vlly = cpl_parameter_get_int(p));
00284         check_nomsg(p = cpl_parameterlist_find(config, "sinfoni.objnod.vurx"));
00285         check_nomsg(vurx = cpl_parameter_get_int(p));
00286         check_nomsg(p = cpl_parameterlist_find(config, "sinfoni.objnod.vury"));
00287         check_nomsg(vury = cpl_parameter_get_int(p));
00288        cknull(cube_tmp = (cpl_imagelist**) cpl_calloc (cfg->nframes, 
00289                                                       sizeof (cpl_imagelist*)),
00290                                      "Could not allocate memory for cube_tmp");
00291 
00292         cknull(cubeobject = (cpl_imagelist**) cpl_calloc (cfg->nframes,
00293                                                   sizeof (cpl_imagelist*)),
00294                                   "Could not allocate memory for cubeobject");
00295         check_nomsg(p=cpl_parameterlist_find(config,
00296                     "sinfoni.sinfo_utl_skycor.rot_cor"));
00297         check_nomsg(sky_cor=cpl_parameter_get_bool(p));
00298 
00299     for ( n = 0 ; n < cfg->nframes ; n++ ) {
00300       if(sky_cor == 1 && pdensity == 2) {
00301         snprintf(file_name,MAX_NAME_SIZE-1,"%s%2.2d%s","out_cube_obj_cor",
00302              n,".fits");
00303 
00304       } else {
00305         snprintf(file_name,MAX_NAME_SIZE-1,"%s%2.2d%s","out_cube_obj",n,
00306              ".fits");
00307       }
00308       check_nomsg(cube_tmp[n] = cpl_imagelist_load(file_name,
00309                                CPL_TYPE_FLOAT,0));
00310       check_nomsg(cubeobject[n] = sinfo_new_cube_getvig(cube_tmp[n],
00311                                 1+vllx,1+vlly,
00312                                 64-vurx,64-vury));
00313           check_nomsg(sinfo_free_imagelist(&cube_tmp[n]));
00314     }
00315         sinfo_free_array_imagelist(&cube_tmp);
00316 
00317     }
00318 
00319     /*
00320         ---------------------------------------------------------------------
00321         ------------------------JITTERING------------------------------------
00322         ---------------------------------------------------------------------
00323     */
00324 
00325     if (cfg->jitterind == 1)
00326     {
00327     sinfo_msg("Jittering...");
00328 
00329         sinfo_msg("Coadded cube size. x: %d y: %d",
00330              cfg->size_x,cfg->size_y);
00331         check_nomsg(jittercube = cpl_imagelist_new()) ;
00332 
00333 
00334     /* 
00335         ---------------------------------------------------------------------
00336         -------------------THOMAS ALGORITHM----------------------------------
00337         ---------------------------------------------------------------------
00338     */
00339        
00340         check_nomsg(p=cpl_parameterlist_find(config,
00341                                              "sinfoni.objnod.scales_sky"));
00342         check_nomsg(scales_sky=cpl_parameter_get_bool(p));
00343         check_nomsg(p=cpl_parameterlist_find(config,"sinfoni.objnod.ks_clip"));
00344         check_nomsg(ks_clip = cpl_parameter_get_bool(p));
00345         check_nomsg(p=cpl_parameterlist_find(config,"sinfoni.objnod.kappa"));
00346         check_nomsg(kappa = cpl_parameter_get_double(p));
00347 
00348 
00349     if(scales_sky == 1) {
00350       sinfo_msg("Subtract spatial sinfo_median to each cube plane");
00351       for(n=0;n<cfg->nframes;n++) {
00352         sinfo_msg("process cube %d\n",n);
00353         sinfo_new_sinfoni_correct_median_it(&(cubeobject[n]));
00354       }
00355     }
00356 
00357 
00358         /* AMO CHECK */
00359 
00360         cknull(maskcube=cpl_imagelist_new(),"could not allocate cube!");
00361         
00362     /* Illumination correction */ 
00363         if(ill_cor != NULL) {
00364       for(n=0;n<cfg->nframes;n++) {
00365         cpl_imagelist_divide_image(cubeobject[n],ill_cor);
00366       }
00367     }
00368         sinfo_free_image(&ill_cor);
00369     
00370     sinfo_msg("Combine jittered cubes");
00371 
00372 
00373     if(ks_clip == 1) {
00374       sinfo_msg("Cube coaddition with kappa-sigma");
00375     }
00376         check_nomsg(onp=cpl_imagelist_get_size(cubeobject[0]));
00377     for(z=0;z<onp;z+=z_stp) {
00378       z_siz=(z_stp < (onp-z)) ? z_stp : (onp-z);
00379       z_min=z;
00380       z_max=z_min+z_siz;
00381       sinfo_msg("Coadding cubes: range [%4.4d,%4.4d) of 0-%d\n",
00382                            z_min,z_max,onp);
00383 
00384       for(j=z_min;j<z_max;j++) {
00385             check_nomsg(j_img=cpl_image_new(cfg->size_x,
00386                                             cfg->size_y,CPL_TYPE_FLOAT));
00387         check_nomsg(cpl_imagelist_set(jittercube,j_img,j));
00388             check_nomsg(m_img = cpl_image_new(cfg->size_x,
00389                                               cfg->size_y,CPL_TYPE_FLOAT));
00390             check_nomsg(cpl_imagelist_set(maskcube,m_img,j));
00391       }
00392       if(ks_clip == 1) {
00393         sinfo_new_combine_jittered_cubes_thomas_range(cubeobject,
00394                             jittercube,
00395                             maskcube,
00396                             cfg->nframes,
00397                             offsetx,offsety,
00398                             times,
00399                             cfg->kernel_type,
00400                             z_min,
00401                             z_max,
00402                             kappa);
00403 
00404       } else {
00405         sinfo_new_combine_jittered_cubes_range(cubeobject, 
00406                                              jittercube,
00407                                              maskcube,
00408                          cfg->nframes,
00409                                              offsetx,
00410                          offsety, 
00411                                              times,
00412                          cfg->kernel_type,
00413                                              z_min,
00414                                              z_max) ;
00415       }
00416     }
00417         sinfo_new_convert_0_to_ZERO_for_cubes(jittercube) ; 
00418  
00419     if (jittercube == NULL)
00420     {
00421         sinfo_msg_error(" could not allocate new data cube!") ;
00422         goto cleanup;
00423     }
00424 
00425         if (maskcube == NULL)
00426     {
00427         sinfo_msg_error(" could not merge the jittered data cubes\n") ;
00428             goto cleanup;
00429     }
00430 
00431     for ( i = 0 ; i <cfg->nframes  ; i++ ) {
00432       sinfo_free_imagelist(&cubeobject[i]);
00433     }
00434     sinfo_free_array_imagelist(&cubeobject);
00435 
00436         newcenter_x = cfg->size_x / 2. + 0.5 ;
00437     newcenter_y = cfg->size_y / 2. + 0.5 ;
00438 
00439         ck0(sinfo_pro_save_ims(jittercube,sof,sof,cfg->outName,
00440                    procatg,NULL,plugin_id,config),
00441         "cannot save cube %s", cfg->outName);
00442 
00443         /* we need to set again the following 3 */
00444         plist=cpl_propertylist_load(file_name,0);
00445         dis=sinfo_pfits_get_cdelt3(plist);
00446         centralLambda=sinfo_pfits_get_crval3(plist);
00447         centralpix=sinfo_pfits_get_crpix3(plist);
00448         sinfo_free_propertylist(&plist);
00449     sinfo_new_set_wcs_cube(jittercube, cfg->outName, 
00450                                centralLambda, dis, 
00451                                centralpix, center_x, center_y);
00452 
00453         cknull(jitter_image = sinfo_new_median_cube(jittercube),
00454                " could not do sinfo_medianCube()");
00455 
00456 
00457         ck0(sinfo_pro_save_ima(jitter_image,sof,sof,cfg->med_cube_name,
00458                    pro_med,NULL,plugin_id,config),
00459         "cannot save ima %s", cfg->outName);
00460 
00461     sinfo_new_set_wcs_image(jitter_image, cfg->med_cube_name,
00462                                 center_x,center_y);
00463 
00464         sinfo_free_image(&jitter_image);
00465 
00466         ck0(sinfo_pro_save_ims(maskcube,sof,sof,cfg->maskname,
00467                    pro_mjit,NULL,plugin_id,config),
00468         "cannot save cube %s", cfg->maskname);
00469 
00470     sinfo_new_set_wcs_cube(maskcube, cfg->maskname, 
00471                                centralLambda, dis, centralpix, 
00472                                center_x, center_y);
00473 
00474         sinfo_free_double(&times) ;
00475         sinfo_free_float(&offsetx) ;
00476         sinfo_free_float(&offsety) ;
00477         sinfo_free_imagelist(&maskcube) ;
00478         sinfo_free_imagelist(&jittercube) ;  
00479  
00480     } /* end of jittering */
00481 
00482  exit:
00483 
00484  /* free memory */
00485     sinfo_objnod_free(&cfg);
00486     sinfo_free_frameset(&stk);
00487     return 0;   
00488 
00489  cleanup:
00490 
00491     sinfo_free_image(&jitter_image);
00492     sinfo_free_imagelist(&jittercube) ;
00493     sinfo_free_imagelist(&maskcube) ;
00494 
00495     if(cfg != NULL) {
00496       if(cube_tmp != NULL) {
00497     for ( n = 0 ; n < cfg->nframes ; n++ ) {
00498       sinfo_free_imagelist(&(cube_tmp[n]));
00499     }
00500         sinfo_free_array_imagelist(&cube_tmp);
00501       }
00502       if(cubeobject != NULL) {
00503     for ( n = 0 ; n < cfg->nframes ; n++ ) {
00504       sinfo_free_imagelist(&(cubeobject[n]));
00505     }
00506         sinfo_free_array_imagelist(&cubeobject);
00507       }
00508 
00509     }
00510     sinfo_free_table(&qclog_tbl);
00511     sinfo_free_image(&im);
00512     sinfo_free_image(&ill_cor);
00513     sinfo_free_float(&offsety);
00514     sinfo_free_float(&offsetx);
00515     sinfo_free_double(&times);
00516     sinfo_objnod_free(&cfg);
00517     sinfo_free_frameset(&stk);
00518     irplib_error_dump(CPL_MSG_ERROR,CPL_MSG_ERROR);
00519     return -1;   
00520 
00521 
00522 
00523 }
00524 

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