sinfo_new_objnod.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_objnod.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_objnod.h"
00036 #include "sinfo_pro_save.h"
00037 #include "sinfo_objnod_ini_by_cpl.h"
00038 #include "sinfo_functions.h"
00039 #include "sinfo_pfits.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 
00064 /*----------------------------------------------------------------------------
00065    Function     : sinfo_new_objnod()
00066    In           : ini_file: file name of according .ini file
00067    Out          : integer (0 if it worked, -1 if it doesn't) 
00068    Job          : this routine carries through the data cube creation of an 
00069                   object science observation using object-sky nodding
00070                   and jittering. This script expects jittered frames that
00071           were already sky-subtracted
00072                   averaged, flatfielded, spectral tilt corrected and 
00073         interleaved if necessary
00074  ---------------------------------------------------------------------------*/
00075 int sinfo_new_objnod (const char* plugin_id,cpl_parameterlist* config, 
00076             cpl_frameset* sof, const char* procatg)
00077 {
00078 
00079     object_config * cfg=NULL ;
00080     cpl_image * im=NULL ;
00081     cpl_image * wavemapim=NULL ;
00082     cpl_image * resampledImage=NULL ;
00083     cpl_image * calim=NULL ;
00084     cpl_image * halospec=NULL ;
00085     cpl_image * sky_im=NULL;
00086     cpl_image* res_flat=NULL;
00087     cpl_image* res_sky=NULL;
00088     cpl_image* flat_im=NULL;
00089     cpl_image* jitter_image=NULL;
00090     cpl_image* eima_avg=NULL;
00091     cpl_image* eima_med=NULL;
00092     cpl_imagelist  * cube=NULL ;
00093     cpl_imagelist  * outcube=NULL ;
00094     cpl_imagelist  * outcube2=NULL ;
00095     cpl_imagelist  ** cubeobject=NULL ;
00096     cpl_imagelist  ** cube_tmp=NULL ;
00097     cpl_imagelist  * jittercube=NULL ;
00098     cpl_imagelist  * maskcube=NULL ;
00099     cpl_imagelist* cflat=NULL;
00100     cpl_imagelist* cflat2=NULL;
00101     cpl_imagelist* csky=NULL;
00102     cpl_imagelist* csky2=NULL;
00103 
00104 
00105     int i=0;
00106     int n=0;
00107     int partind = 0 ;
00108     int centralpix=0 ;
00109     int z_siz=0;
00110     int z_min=0;
00111     int z_max=0;
00112     int z=0;
00113     int z_stp=100;
00114     int scales_sky=0;
00115     int ks_clip=0;
00116     double kappa=2.0;
00117 
00118     float ref_offx=0;
00119     float ref_offy=0;
00120     float mi=0 ;
00121     float ma=0 ;
00122     float fcol=0 ;
00123     float center_x=0;
00124     float newcenter_x=0 ;
00125     float center_y=0;
00126     float newcenter_y=0;
00127     float cd1_1=0;
00128     float cd1_2=0;
00129     float cd2_1=0;
00130     float cd2_2=0;
00131     float pixelscale=0;
00132     float angle=0;
00133     float radangle=0;
00134     double exptime=0;
00135 
00136     float *  correct_dist=NULL ;
00137     float * distances=NULL ;
00138     double * times=NULL ;
00139     float * offsetx=NULL;
00140     float * offsety=NULL;
00141     float ** slit_edges=NULL ;
00142 
00143     double dis=0;
00144     double centralLambda=0;
00145 
00146     char name_jitter[MAX_NAME_SIZE] ;
00147     char pro_mjit[MAX_NAME_SIZE];
00148     char pro_obs[MAX_NAME_SIZE];
00149     char pro_med[MAX_NAME_SIZE];
00150 
00151 
00152     char * name=NULL ;
00153     char * partname=NULL;
00154     char * partname2=NULL;
00155     char file_name[MAX_NAME_SIZE];
00156     int vllx=0;
00157     int vlly=0;
00158     int vurx=0;
00159     int vury=0;
00160 
00161     int onp=0;
00162     int j=0;
00163     cpl_image* j_img=NULL;
00164     cpl_image* m_img=NULL;
00165     cpl_table* qclog_tbl=NULL;
00166     cpl_image* ill_cor=NULL;
00167     cpl_frame* frame=NULL;
00168     cpl_frameset* stk=NULL;
00169     cpl_parameter* p=NULL;
00170     cpl_propertylist* plist=NULL;
00171 
00172      if (strcmp(procatg,PRO_COADD_STD) == 0) {
00173     strcpy(pro_mjit,PRO_MASK_COADD_STD);
00174     strcpy(pro_obs,PRO_OBS_STD);
00175     strcpy(pro_med,PRO_MED_COADD_STD);
00176 
00177     } else if (strcmp(procatg,PRO_COADD_PSF) == 0) {
00178     strcpy(pro_mjit,PRO_MASK_COADD_PSF);
00179     strcpy(pro_obs,PRO_OBS_PSF);
00180     strcpy(pro_med,PRO_MED_COADD_PSF);
00181     } else {
00182     strcpy(pro_mjit,PRO_MASK_COADD_OBJ);
00183     strcpy(pro_obs,PRO_OBS_OBJ);
00184     strcpy(pro_med,PRO_MED_COADD_OBJ);
00185     }
00186 
00187 
00188     /*----parse input data and parameters to set cube_config cfg---*/
00189     check_nomsg(stk = cpl_frameset_new());
00190 
00191     cknull(cfg = sinfo_parse_cpl_input_objnod(config,sof,&stk),
00192        "Error setting parameter configuration");
00193  
00194     ck0(sinfo_check_input_data(cfg),"error checking input");
00195 
00196     if ( cfg->jitterind == 1 )
00197     {
00198         cknull(times = (double*) cpl_calloc (cfg->nframes, sizeof (double)), 
00199            " could not allocate memory!") ;
00200  
00201         cknull(offsetx = (float*) cpl_calloc (cfg->nframes, sizeof(float)),
00202            " could not allocate memory!") ;
00203 
00204         cknull(offsety = (float*) cpl_calloc (cfg->nframes, sizeof(float)),
00205            " could not allocate memory!") ;
00206     }
00207   
00208     if (cfg->jitterind == 0)
00209     {
00210         if ( NULL != (partname = strtok(cfg->outName, "."))) 
00211         {
00212             partname2 = strtok (NULL, ".") ;
00213             partind = 1 ;
00214         }
00215     }
00216 
00217 
00218     ck0(sinfo_auto_size_cube2(cfg,&ref_offx,&ref_offy),"Error resizing cube");
00219   
00220 
00221 
00222     check_nomsg(p=cpl_parameterlist_find(config,"sinfoni.objnod.fcol"));
00223     check_nomsg(fcol=cpl_parameter_get_double(p));
00224    
00225     if(NULL != cpl_frameset_find(sof,PRO_ILL_COR)) {
00226      sinfo_msg("Illumination correction cube is provided");
00227      frame = cpl_frameset_find(sof,PRO_ILL_COR);
00228      ill_cor=cpl_image_load(cpl_frame_get_filename(frame),CPL_TYPE_FLOAT,0,0);
00229     } else {
00230      sinfo_msg("Illumination correction cube not provided");
00231      cpl_error_reset();
00232     }
00233 
00234      for ( n = 0 ; n < cfg->nframes ; n++ )
00235      {
00236 
00237         sinfo_msg("Read FITS information");
00238         name = cfg->framelist[n] ;    
00239         if (n == 0)
00240     {
00241         strcpy (name_jitter, name) ;
00242     }   
00243         if( sinfo_is_fits_file(name) != 1) {
00244           sinfo_msg_error("Input file %s is not FITS",name);
00245           goto cleanup;
00246     }
00247     
00248     /* get some header values and compute the CD-sinfo_matrix */
00249         plist=cpl_propertylist_load(name,0);
00250     pixelscale = sinfo_pfits_get_pixscale(plist) /2;
00251     angle = sinfo_pfits_get_posangle(plist) ;
00252         /* in PUPIL data there is not posangle info: we reset the error */
00253         if(cpl_error_get_code() != CPL_ERROR_NONE) {
00254       cpl_error_reset();
00255         }
00256     sinfo_free_propertylist(&plist);
00257     radangle = angle * PI_NUMB / 180. ;
00258     cd1_1 = cos(radangle) ;
00259     cd1_2 = sin(radangle) ;
00260     cd2_1 = -sin(radangle) ;
00261     cd2_2 = cos(radangle) ;
00262     
00263     sinfo_msg("frame no.: %d, name: %s\n", n, name) ;
00264     cknull(im = cpl_image_load(name,CPL_TYPE_FLOAT,0,0),
00265            " could not load frame %s!",name) ;
00266 
00267     if (cfg->jitterind == 1)
00268     {
00269       exptime = sinfo_pfits_get_ditndit(name) ;
00270 
00271       if (exptime == FLAG)
00272         {
00273           sinfo_msg_error("could not read fits header keyword exptime!");
00274           goto cleanup;
00275         }
00276       times[n] = exptime ;
00277           ck0(sinfo_new_assign_offset(n,name,offsetx,offsety,
00278                                       ref_offx,ref_offy),
00279                                       "Error assigning offsets");
00280 
00281     }
00282 
00283         /*
00284          *--------------------------------------------------------------    
00285          *---------------------RESAMPLING-------------------------------
00286      *--------------------------------------------------------------
00287          */
00288         sinfo_msg("Resampling object");
00289     cknull(wavemapim = cpl_image_load(cfg->wavemap,CPL_TYPE_FLOAT,0,0),
00290            "could not load wavemap");
00291 
00292         cknull(resampledImage = sinfo_new_defined_resampling( im, 
00293                                   wavemapim, 
00294                                   cfg->ncoeffs,
00295                                   &cfg->nrows,
00296                                   &dis,
00297                                   &mi,
00298                                   &ma,
00299                                   &centralLambda,
00300                                   &centralpix),
00301            " sinfo_definedResampling() failed" ) ;
00302 
00303  
00304     if(n ==0) {
00305           if(strcmp(cfg->mflat_dist,"not_found") != 0) {
00306         sinfo_msg("Resampling master flat");
00307         cknull(flat_im=cpl_image_load(cfg->mflat_dist,CPL_TYPE_FLOAT,0,0),
00308                    "Distorted master flat field not found\n"
00309                    "You may have set --stack-flat_ind=FALSE\n"
00310                    "Flat field resampling skipped");
00311             cknull(res_flat = sinfo_new_defined_resampling(flat_im, 
00312                                wavemapim,
00313                                cfg->ncoeffs,
00314                                &cfg->nrows,
00315                                &dis,
00316                                &mi,
00317                                &ma,
00318                                &centralLambda,
00319                                &centralpix),
00320            " sinfo_definedResampling() failed" ) ;
00321 
00322             sinfo_free_image(&flat_im) ;
00323       }
00324      
00325       if(strcmp(cfg->sky_dist,"no_sky")!=0) {
00326         sinfo_msg("Resampling sky");
00327         check_nomsg(sky_im=cpl_image_load(cfg->sky_dist,
00328                                               CPL_TYPE_FLOAT,0,0));
00329         cknull(res_sky = sinfo_new_defined_resampling(sky_im,
00330                               wavemapim,
00331                               cfg->ncoeffs,
00332                               &cfg->nrows,
00333                               &dis,
00334                               &mi,
00335                               &ma,
00336                               &centralLambda,
00337                               &centralpix),
00338            " sinfo_definedResampling() failed" );
00339 
00340         sinfo_free_image(&sky_im) ;
00341 
00342        
00343       }
00344     }
00345 
00346         sinfo_msg ("dispersion %f\n", dis) ;
00347         sinfo_msg ("lambda min %f max %f cent %f\n", mi,ma,centralLambda ) ;
00348         sinfo_msg ("central pixel %d\n", centralpix) ; 
00349 
00350         sinfo_free_image(&im) ;
00351         sinfo_free_image(&wavemapim) ;
00352 
00353     /*
00354          *-------------------------------------------------------------------
00355          *----------------Calibration----------------------------------------
00356          *-------------------------------------------------------------------
00357     */
00358         /*----Multiply with calibrated halogen lamp spectrum----*/ 
00359         if (cfg->halocorrectInd == 1)
00360         {     
00361       sinfo_msg("Calibration");
00362       check_nomsg(halospec = cpl_image_load(cfg->halospectrum,
00363                                                 CPL_TYPE_FLOAT,0,0)) ;
00364 
00365       cknull(calim = sinfo_new_multiply_image_with_spectrum(resampledImage,
00366                                                                 halospec),
00367                         " sinfo_new_multiply_image_with_spectrum() failed" ) ;
00368 
00369       sinfo_free_image(&halospec) ;
00370       sinfo_free_image(&resampledImage) ;
00371       resampledImage = cpl_image_duplicate(calim) ;
00372       sinfo_free_image(&calim);
00373         }
00374     
00375         /*
00376          *-------------------------------------------------------------------
00377          *------------------CUBECREATION-------------------------------------
00378          *-------------------------------------------------------------------
00379      */
00380     sinfo_msg("Cube creation");
00381     /*---select north-south-test or fitting of slitlet edges--*/
00382     if (cfg->northsouthInd == 0) {
00383       sinfo_msg("cfg->northsouthInd == 0");
00384       cknull(slit_edges = sinfo_read_slitlets_edges(cfg),
00385          "error reading slitlets edges");
00386     } else {
00387       sinfo_msg("cfg->northsouthInd != 0");
00388       cknull(distances = sinfo_read_distances(cfg),
00389          "error reading distances");
00390     }
00391  
00392     cknull(correct_dist = (float*) cpl_calloc(cfg->nslits, sizeof (float)),
00393            " could not allocate memory!") ;
00394 
00395     sinfo_msg("Create cube object");
00396     if (cfg->northsouthInd ==0 ) {
00397 
00398         cknull(cube = sinfo_new_make_cube_spi(resampledImage,
00399                                                   slit_edges,
00400                                                   correct_dist),
00401                                  " could not construct data cube!") ;
00402 
00403     }  else {
00404        cknull(cube = sinfo_new_make_cube_dist(resampledImage,
00405                                                   fcol,
00406                                                   distances,
00407                                                   correct_dist),
00408                                          " could not construct a data cube!") ;
00409     }
00410     sinfo_free_image(&resampledImage);
00411 
00412     if(n==0) {
00413       if(strcmp(cfg->mflat_dist,"not_found")!=0) {
00414         sinfo_msg("Create cube master flat");
00415         if (cfg->northsouthInd ==0 ) {
00416           cknull(cflat=sinfo_new_make_cube_spi(res_flat,
00417                                                    slit_edges,
00418                                                    correct_dist),
00419                                 " could not construct data cube!") ;
00420         }  else {
00421           cknull(cflat = sinfo_new_make_cube_dist(res_flat,
00422                                                       fcol,
00423                                                       distances,
00424                                                       correct_dist),
00425                                         " could not construct a data cube!") ;
00426         }
00427         sinfo_free_image(&res_flat);
00428       }
00429       if(strcmp(cfg->sky_dist,"no_sky")!=0) {
00430 
00431         sinfo_msg("Create cube sky");
00432         if (cfg->northsouthInd ==0 ) {
00433           cknull(csky = sinfo_new_make_cube_spi(res_sky,
00434                                                     slit_edges,
00435                                                     correct_dist),
00436                                            " could not construct data cube!") ;
00437         }  else {
00438           cknull(csky = sinfo_new_make_cube_dist(res_sky,
00439                                                      fcol,
00440                                                      distances,
00441                                                      correct_dist),
00442                                          " could not construct a data cube!") ;
00443         }
00444         sinfo_free_image(&res_sky);
00445       }
00446     }
00447 
00448 
00449         if (cfg->northsouthInd ==0 )
00450      {
00451        sinfo_new_destroy_2Dfloatarray(&slit_edges,cfg->nslits);
00452      }
00453     else
00454       {
00455             sinfo_new_destroy_array(&distances);
00456       }
00457 
00458         /*
00459          *--------------------------------------------------------------------
00460          *------------------------FINETUNING----------------------------------
00461          *--------------------------------------------------------------------
00462          * shift the rows of the reconstructed images of the data cube to the 
00463          * correct sub pixel position select the shift method: polynomial 
00464          * interpolation, FFT or cubic spline interpolation
00465          *--------------------------------------------------------------------
00466          */
00467  
00468     if(n==0) {
00469       if(strcmp(cfg->sky_dist,"no_sky")!=0) {
00470          cknull(csky2=sinfo_new_fine_tune(csky,
00471                                              correct_dist,
00472                                              cfg->method,
00473                                              cfg->order,
00474                           cfg->nslits),
00475             " could not fine tune the data cube") ;
00476 
00477          sinfo_free_imagelist(&csky);
00478          sinfo_msg("Stretch output cube along Y direction");
00479  
00480          cknull(csky = sinfo_new_bin_cube(csky2,1,2,0,63,0,63),
00481             "error rebinning sky cube");
00482          sinfo_free_imagelist(&csky2);
00483 
00484        
00485      
00486          ck0(sinfo_pro_save_ims(csky,sof,sof,"out_sky_cube.fits",
00487                     PRO_OBS_SKY,NULL,plugin_id,config),
00488          "cannot dump cube %s", "out_sky_cube.fits");
00489  
00490          cknull(eima_med=sinfo_new_median_cube(csky),
00491                 "Creating an average image");
00492          check_nomsg(center_x = cpl_image_get_size_x(eima_med)/ 2. + 0.5);
00493          check_nomsg(center_y = cpl_image_get_size_y(eima_med)/ 2. + 0.5);
00494 
00495          sinfo_new_set_wcs_cube(csky, "out_sky_cube.fits", centralLambda, 
00496                           dis, centralpix, center_x, center_y);
00497 
00498          sinfo_free_imagelist(&csky) ;
00499 
00500          ck0(sinfo_pro_save_ima(eima_med,sof,sof,"out_sky_med.fits",
00501                     PRO_SKY_MED,NULL,plugin_id,config),
00502          "cannot save ima %s", "out_sky_med.fits");
00503 
00504          sinfo_new_set_wcs_image(eima_med,"out_sky_med.fits", 
00505                                      center_x, center_y);
00506          sinfo_free_image(&eima_med);
00507       }
00508 
00509 
00510       if(strcmp(cfg->mflat_dist,"not_found")!=0) {
00511 
00512         cknull(cflat2=sinfo_new_fine_tune(cflat,correct_dist,
00513                                               cfg->method,cfg->order,
00514                                               cfg->nslits),
00515                                        " could not fine tune the data cube") ;
00516 
00517         sinfo_free_imagelist(&cflat);
00518         sinfo_msg("Stretch output cube along Y direction");
00519  
00520         cknull(cflat = sinfo_new_bin_cube(cflat2,1,2,0,63,0,63),
00521            "Error binning flat cube");
00522         sinfo_free_imagelist(&cflat2);
00523 
00524         ck0(sinfo_pro_save_ims(cflat,sof,sof,"out_mflat_cube.fits",
00525                    "MFLAT_CUBE",NULL,plugin_id,config),
00526         "cannot save cube %s", "out_mflat_cube.fits");
00527 
00528         cknull(eima_avg=sinfo_new_average_cube_to_image(cflat),
00529                "Creating an average image");
00530 
00531         ck0(sinfo_pro_save_ima(eima_avg,sof,sof,"out_mflat_avg.fits",
00532                    "MFLAT_AVG",NULL,plugin_id,config),
00533         "cannot save ima %s", "out_mflat_avg.fits");
00534 
00535         sinfo_free_image(&eima_avg);
00536 
00537         cknull(eima_med=sinfo_new_median_cube(cflat),
00538            "Error computing median on cube flat");
00539 
00540         ck0(sinfo_pro_save_ima(eima_med,sof,sof,"out_mflat_med.fits",
00541                    "MFLAT_MED",NULL,plugin_id,config),
00542         "cannot save ima %s", "out_mflat_med.fits");
00543 
00544         sinfo_free_imagelist(&cflat); 
00545         sinfo_free_image(&eima_med);
00546       }
00547     }
00548        
00549 
00550         cknull(outcube2=sinfo_new_fine_tune(cube,
00551                                             correct_dist,
00552                                             cfg->method,
00553                                             cfg->order,
00554                                             cfg->nslits),
00555                                      " could not fine tune the data cube") ;
00556 
00557         sinfo_msg("Stretch output cube along Y direction");
00558         cknull(outcube = sinfo_new_bin_cube(outcube2,1,2,0,63,0,63),
00559            "Error binning cube");
00560         sinfo_free_imagelist(&cube);
00561         cknull_nomsg(qclog_tbl=sinfo_qclog_init());
00562         sinfo_get_pupil_shift(outcube,n,&qclog_tbl);
00563         snprintf(file_name,MAX_NAME_SIZE-1,"%s%2.2d%s","out_cube_obj",n,".fits");
00564         ck0(sinfo_pro_save_ims(outcube,sof,sof,file_name,
00565                    pro_obs,qclog_tbl,plugin_id,config),
00566         "cannot save cube %s", file_name);
00567 
00568 
00569         sinfo_free_table(&qclog_tbl);
00570     check_nomsg(center_x = cpl_image_get_size_x(
00571                                cpl_imagelist_get(outcube,0))/2.+0.5) ;
00572     check_nomsg(center_y = cpl_image_get_size_y(
00573                                cpl_imagelist_get(outcube,0))/2.+0.5 );
00574    
00575 
00576     sinfo_new_set_wcs_cube(outcube, file_name, centralLambda, dis, 
00577                      centralpix, center_x, center_y);
00578 
00579     /* free memory */
00580         /* to prevent error message comment next line */
00581         sinfo_free_imagelist(&outcube2);
00582         sinfo_free_imagelist(&outcube) ;
00583     sinfo_free_float(&correct_dist) ;
00584 
00585 
00586     } /* end loop over n (nframes) */
00587  
00588     /* leak free */
00589     if(cfg->jitterind == 0) {
00590       goto exit;
00591     }   
00592      
00593     /* Here in case of autojitter we estimate the sky */
00594     if( (cfg->size_x*cfg->size_y*cfg->nframes) > (130*130*21) ) {
00595       sinfo_msg_warning("Coadd cube size:%d,%d. N frames: %d",
00596                            cfg->size_x,cfg->size_y,cfg->nframes);
00597       sinfo_msg_warning("Max allowed should be such that "
00598                         "sixeX*sixeY*Nframes < 130*130*21");
00599       goto exit;
00600     } 
00601  
00602     if ( cfg->jitterind == 1 )
00603     {  
00604         check_nomsg(p = cpl_parameterlist_find(config, "sinfoni.objnod.vllx"));
00605         check_nomsg(vllx = cpl_parameter_get_int(p));
00606         check_nomsg(p = cpl_parameterlist_find(config, "sinfoni.objnod.vlly"));
00607         check_nomsg(vlly = cpl_parameter_get_int(p));
00608         check_nomsg(p = cpl_parameterlist_find(config, "sinfoni.objnod.vurx"));
00609         check_nomsg(vurx = cpl_parameter_get_int(p));
00610         check_nomsg(p = cpl_parameterlist_find(config, "sinfoni.objnod.vury"));
00611         check_nomsg(vury = cpl_parameter_get_int(p));
00612  
00613         cknull(cube_tmp =(cpl_imagelist**) cpl_calloc(cfg->nframes, 
00614                                                       sizeof (cpl_imagelist*)),
00615                                      "Could not allocate memory for cube_tmp");
00616         cknull(cubeobject =(cpl_imagelist**) cpl_calloc(cfg->nframes, 
00617                                                         sizeof(cpl_imagelist*)),
00618                                     "Could not allocate memory for cubeobject");
00619 
00620     for ( n = 0 ; n < cfg->nframes ; n++ ) {
00621           snprintf(file_name,MAX_NAME_SIZE-1,"%s%2.2d%s","out_cube_obj",n,".fits");
00622       check_nomsg(cube_tmp[n] = cpl_imagelist_load(file_name,
00623                                CPL_TYPE_FLOAT,0));
00624       check_nomsg(cubeobject[n] = sinfo_new_cube_getvig(cube_tmp[n],
00625                                 1+vllx,1+vlly,
00626                                 64-vurx,64-vury));
00627           check_nomsg(sinfo_free_imagelist(&cube_tmp[n]));
00628     }
00629         sinfo_free_array_imagelist(&cube_tmp);
00630 
00631     }
00632 
00633     /*
00634         ---------------------------------------------------------------------
00635         ------------------------JITTERING------------------------------------
00636         ---------------------------------------------------------------------
00637     */
00638 
00639     if (cfg->jitterind == 1)
00640     {
00641     sinfo_msg("Jittering...");
00642 
00643         sinfo_msg("Coadded cube size. x: %d y: %d",
00644              cfg->size_x,cfg->size_y);
00645         check_nomsg(jittercube = cpl_imagelist_new()) ;
00646 
00647 
00648     /* 
00649         ---------------------------------------------------------------------
00650         -------------------THOMAS ALGORITHM----------------------------------
00651         ---------------------------------------------------------------------
00652     */
00653        
00654         check_nomsg(p=cpl_parameterlist_find(config,
00655                                              "sinfoni.objnod.scales_sky"));
00656         check_nomsg(scales_sky=cpl_parameter_get_bool(p));
00657         check_nomsg(p=cpl_parameterlist_find(config,"sinfoni.objnod.ks_clip"));
00658         check_nomsg(ks_clip = cpl_parameter_get_bool(p));
00659         check_nomsg(p=cpl_parameterlist_find(config,"sinfoni.objnod.kappa"));
00660         check_nomsg(kappa = cpl_parameter_get_double(p));
00661 
00662 
00663     if(scales_sky == 1) {
00664       sinfo_msg("Subtract spatial sinfo_median to each cube plane");
00665       for(n=0;n<cfg->nframes;n++) {
00666         sinfo_msg("process cube %d\n",n);
00667         sinfo_new_sinfoni_correct_median_it(&(cubeobject[n]));
00668       }
00669     }
00670 
00671 
00672         /* AMO CHECK */
00673 
00674         cknull(maskcube=cpl_imagelist_new(),"could not allocate cube!");
00675         
00676     /* Illumination correction */ 
00677         if(ill_cor != NULL) {
00678       for(n=0;n<cfg->nframes;n++) {
00679         sinfo_msg("Illumination correction is applied");
00680         cpl_imagelist_divide_image(cubeobject[n],ill_cor);
00681       }
00682     }
00683         sinfo_free_image(&ill_cor);
00684     
00685     sinfo_msg("Combine jittered cubes");
00686 
00687 
00688     if(ks_clip == 1) {
00689       sinfo_msg("Cube coaddition with kappa-sigma");
00690     }
00691         check_nomsg(onp=cpl_imagelist_get_size(cubeobject[0]));
00692     for(z=0;z<onp;z+=z_stp) {
00693       z_siz=(z_stp < (onp-z)) ? z_stp : (onp-z);
00694       z_min=z;
00695       z_max=z_min+z_siz;
00696       sinfo_msg("Coadding cubes: range [%4.4d,%4.4d) of 0-%d\n",
00697                            z_min,z_max,onp);
00698 
00699       for(j=z_min;j<z_max;j++) {
00700             check_nomsg(j_img=cpl_image_new(cfg->size_x,
00701                                             cfg->size_y,CPL_TYPE_FLOAT));
00702         check_nomsg(cpl_imagelist_set(jittercube,j_img,j));
00703             check_nomsg(m_img = cpl_image_new(cfg->size_x,
00704                                               cfg->size_y,CPL_TYPE_FLOAT));
00705             check_nomsg(cpl_imagelist_set(maskcube,m_img,j));
00706       }
00707       if(ks_clip == 1) {
00708         sinfo_new_combine_jittered_cubes_thomas_range(cubeobject,
00709                             jittercube,
00710                             maskcube,
00711                             cfg->nframes,
00712                             offsetx,offsety,
00713                             times,
00714                             cfg->kernel_type,
00715                             z_min,
00716                             z_max,
00717                             kappa);
00718 
00719       } else {
00720         sinfo_new_combine_jittered_cubes_range(cubeobject, 
00721                                              jittercube,
00722                                              maskcube,
00723                          cfg->nframes,
00724                                              offsetx,
00725                          offsety, 
00726                                              times,
00727                          cfg->kernel_type,
00728                                              z_min,
00729                                              z_max) ;
00730       }
00731     }
00732         sinfo_new_convert_0_to_ZERO_for_cubes(jittercube) ; 
00733 
00734     if (jittercube == NULL)
00735     {
00736         sinfo_msg_error(" could not allocate new data cube!") ;
00737         goto cleanup;
00738     }
00739 
00740         if (maskcube == NULL)
00741     {
00742         sinfo_msg_error(" could not merge the jittered data cubes\n") ;
00743             goto cleanup;
00744     }
00745 
00746     for ( i = 0 ; i <cfg->nframes  ; i++ ) {
00747       sinfo_free_imagelist(&cubeobject[i]);
00748     }
00749     sinfo_free_array_imagelist(&cubeobject);
00750 
00751         newcenter_x = cfg->size_x / 2. + 0.5 ;
00752     newcenter_y = cfg->size_y / 2. + 0.5 ;
00753 
00754         ck0(sinfo_pro_save_ims(jittercube,sof,sof,cfg->outName,
00755                    procatg,NULL,plugin_id,config),
00756         "cannot save cube %s", cfg->outName);
00757 
00758     sinfo_new_set_wcs_cube(jittercube, cfg->outName, centralLambda, 
00759                                dis, centralpix, center_x, center_y);
00760 
00761         cknull(jitter_image = sinfo_new_median_cube(jittercube),
00762                " could not do sinfo_medianCube()");
00763 
00764 
00765         ck0(sinfo_pro_save_ima(jitter_image,sof,sof,cfg->med_cube_name,
00766                    pro_med,NULL,plugin_id,config),
00767         "cannot save ima %s", cfg->outName);
00768 
00769     sinfo_new_set_wcs_image(jitter_image, cfg->med_cube_name,
00770                                 center_x,center_y);
00771 
00772         sinfo_free_image(&jitter_image);
00773 
00774         ck0(sinfo_pro_save_ims(maskcube,sof,sof,cfg->maskname,
00775                    pro_mjit,NULL,plugin_id,config),
00776         "cannot save cube %s", cfg->maskname);
00777 
00778     sinfo_new_set_wcs_cube(maskcube, cfg->maskname, centralLambda, 
00779                                dis, centralpix, center_x, center_y);
00780 
00781         sinfo_free_double(&times) ;
00782         sinfo_free_float(&offsetx) ;
00783         sinfo_free_float(&offsety) ;
00784         sinfo_free_imagelist(&maskcube) ;
00785         sinfo_free_imagelist(&jittercube) ;         
00786 
00787     } /* end of jittering */
00788 
00789  exit:
00790 
00791  /* free memory */
00792     sinfo_objnod_free(&cfg);
00793     sinfo_free_frameset(&stk);
00794     return 0;   
00795 
00796  cleanup:
00797     sinfo_free_propertylist(&plist);
00798     sinfo_free_image(&jitter_image);
00799     sinfo_free_imagelist(&jittercube) ;
00800     sinfo_free_imagelist(&maskcube) ;
00801 
00802     if(cfg != NULL) {
00803       if(cube_tmp != NULL) {
00804     for ( n = 0 ; n < cfg->nframes ; n++ ) {
00805       sinfo_free_imagelist(&(cube_tmp[n]));
00806     }
00807         sinfo_free_array_imagelist(&cube_tmp);
00808       }
00809       if(cubeobject != NULL) {
00810     for ( n = 0 ; n < cfg->nframes ; n++ ) {
00811       sinfo_free_imagelist(&(cubeobject[n]));
00812     }
00813         sinfo_free_array_imagelist(&cubeobject);
00814       }
00815 
00816     }
00817 
00818     sinfo_free_imagelist(&outcube2) ;
00819     sinfo_free_imagelist(&outcube) ;
00820     sinfo_free_table(&qclog_tbl);
00821     sinfo_free_image(&eima_avg);
00822     sinfo_free_image(&eima_med);
00823     sinfo_free_imagelist(&cflat) ;
00824     sinfo_free_imagelist(&cflat2) ;
00825     sinfo_free_imagelist(&cube) ;
00826     sinfo_free_imagelist(&csky) ;
00827     sinfo_free_imagelist(&csky2) ;
00828 
00829     if(cfg!=NULL) {
00830       if (cfg->northsouthInd ==0 ) {
00831     if(slit_edges != NULL) {
00832        sinfo_new_destroy_2Dfloatarray(&slit_edges,cfg->nslits);
00833     }
00834       } else {
00835     if (distances != NULL ) {
00836             sinfo_new_destroy_array(&distances);
00837     }
00838       }
00839     }
00840 
00841     sinfo_free_float(&correct_dist);
00842     sinfo_free_image(&res_flat);
00843     sinfo_free_image(&res_sky);
00844     sinfo_free_image(&calim);
00845     sinfo_free_image(&halospec) ;
00846     sinfo_free_image(&sky_im) ;
00847     sinfo_free_image(&resampledImage);
00848     sinfo_free_image(&flat_im) ;
00849     sinfo_free_image(&wavemapim);
00850     sinfo_free_image(&im);
00851     sinfo_free_image(&ill_cor);
00852     sinfo_free_float(&offsety);
00853     sinfo_free_float(&offsetx);
00854     sinfo_free_double(&times);
00855     sinfo_objnod_free(&cfg);
00856     sinfo_free_frameset(&stk);
00857 
00858     return -1;   
00859 
00860 
00861 
00862 }

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