00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifdef HAVE_CONFIG_H
00029 # include <config.h>
00030 #endif
00031
00032
00033
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
00048
00049 #define PI_NUMB (3.1415926535897932384626433832795)
00050
00051
00052
00053
00054
00055
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
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
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
00224 plist=cpl_propertylist_load(name,0);
00225 pixelscale = sinfo_pfits_get_pixscale(plist) /2;
00226 angle = sinfo_pfits_get_posangle(plist) ;
00227
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 }
00262
00263
00264
00265 if(cfg->jitterind == 0) {
00266 goto exit;
00267 }
00268
00269
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
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
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
00359
00360 cknull(maskcube=cpl_imagelist_new(),"could not allocate cube!");
00361
00362
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
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(×) ;
00475 sinfo_free_float(&offsetx) ;
00476 sinfo_free_float(&offsety) ;
00477 sinfo_free_imagelist(&maskcube) ;
00478 sinfo_free_imagelist(&jittercube) ;
00479
00480 }
00481
00482 exit:
00483
00484
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(×);
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