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
00029
00030
00031 #ifdef HAVE_CONFIG_H
00032 #include <config.h>
00033 #endif
00034
00035
00036
00037
00038
00039
00040 #include <strings.h>
00041 #include <string.h>
00042 #include <stdio.h>
00043 #include <math.h>
00044 #include <libgen.h>
00045
00046
00047
00048 #include <cpl.h>
00049
00050
00051 #include <irplib_utils.h>
00052
00053
00054 #include <sinfo_pro_types.h>
00055 #include <sinfo_product_config.h>
00056 #include <sinfo_prepare_stacked_frames_config.h>
00057 #include <sinfo_objnod_config.h>
00058 #include <sinfo_new_objnod.h>
00059 #include <sinfo_new_prepare_stacked_frames.h>
00060 #include <sinfo_pro_save.h>
00061 #include <sinfo_raw_types.h>
00062 #include <sinfo_functions.h>
00063 #include <sinfo_tpl_utils.h>
00064 #include <sinfo_tpl_dfs.h>
00065 #include <sinfo_hidden.h>
00066 #include <sinfo_globals.h>
00067 #include <sinfo_rec_utils.h>
00068
00069
00070
00071
00072
00073 static int sinfo_utl_illumcorr_create(cpl_plugin *plugin);
00074 static int sinfo_utl_illumcorr_exec(cpl_plugin *plugin);
00075 static int sinfo_utl_illumcorr_destroy(cpl_plugin *plugin);
00076 static int sinfo_utl_illumcorr(cpl_parameterlist *config, cpl_frameset *set);
00077
00078
00079 #define SINFO_DOUBLE_SWAP(a,b) { register double t=(a);(a)=(b);(b)=t; }
00080
00081 static cpl_error_code
00082 sinfo_tools_sort_double(
00083 double * pix_arr,
00084 int n);
00085
00086 static cpl_frame*
00087 sinfo_get_dummy_object(cpl_frameset* obj_set);
00088
00089 static void
00090 sinfo_illumcorr_config_add (cpl_parameterlist *list);
00091
00092 static int
00093 create_illumcorr (const char* plugin_id,
00094 cpl_parameterlist *cpl_cfg,
00095 cpl_frameset* sof,
00096 const char *name_i);
00097 static int
00098 sinfo_illumcorr_create_bins (cpl_imagelist *sky,
00099 int llx, int lly, int urx, int ury,
00100 int spec_bin,
00101 double min_flux,
00102 int ** start,
00103 int ** end,
00104 int z1, int z2);
00105
00106 static int
00107 sinfo_juha_function1d_natural_spline(double *, double *, int, double *,
00108 double *, int);
00109
00110 static int
00111 sinfo_function1d_search_value(double *, int, double, int *) ;
00112
00113 static cpl_vector *
00114 sinfo_vector_filter_median_create(const cpl_vector * v, int hw);
00115
00116 static cpl_vector *
00117 sinfo_juha_vector_filter_median_create(const cpl_vector * v, int hw);
00118
00119 static double
00120 sinfo_image_get_median_window (const cpl_image *image,
00121 int llx, int lly, int urx, int ury);
00122
00123
00124
00125
00126
00127
00128 static char sinfo_utl_illumcorr_description1[] =
00129 "This recipe calculates illumination correction based on sky emission.\n"
00130 "The input files are sky (or object) frames tagged\n"
00131 " SKY_NODDING (OBJECT_NODDING)\n"
00132 "Master calibration frames:\n";
00133
00134
00135 static char sinfo_utl_illumcorr_description2[] =
00136 "A corresponding (DIT) dark frame (tag=MASTER_DARK)"
00137 "A corresponding (band,preoptics) wavelength map image (tag=WAVE_MAP)\n"
00138 "A corresponding (band,preoptics) master flat field (tag=MASTER_FLAT_LAMP)\n"
00139 "A corresponding (band,preoptics) master bad pixel map (tag=MASTER_BP_MAP)\n"
00140 "A corresponding (band,preoptics) slitlets position frame (tag=SLIT_POS)\n"
00141 "A corresponding (band) distortion table (tag=DISTORTION)\n"
00142 "A corresponding (band) slitlet distance table (tag=SLITLETS_DISTANCE)\n";
00143
00144
00145 static char sinfo_utl_illumcorr_description3[] =
00146 "The output is a cube resulting from the analysis of sky emission\n";
00147
00148
00149 static char sinfo_utl_illumcorr_description4[] =
00150 "Information on relevant parameters can be found with\n"
00151 "esorex --params sinfo_utl_illumcorr\n"
00152 "esorex --help sinfo_utl_illumcorr\n"
00153 "\n";
00154
00155 static char sinfo_utl_illumcorr_description[1300];
00156
00157
00158
00159
00160
00161
00166
00168 static int
00169 sinfo_utl_illumcorr_create(cpl_plugin *plugin)
00170 {
00171
00172
00173
00174
00175
00176
00177 cpl_recipe *recipe = (cpl_recipe *)plugin;
00178 recipe->parameters = cpl_parameterlist_new();
00179 if(recipe->parameters == NULL) {
00180 return 1;
00181 }
00182 cpl_error_reset();
00183 irplib_reset();
00184
00185
00186
00187
00188 sinfo_product_config_add (recipe->parameters);
00189 sinfo_prepare_stacked_frames_config_add(recipe->parameters);
00190 sinfo_objnod_config_add(recipe->parameters);
00191 sinfo_illumcorr_config_add (recipe->parameters);
00192
00193 return 0;
00194
00195 }
00196
00197 static int
00198 sinfo_utl_illumcorr_exec(cpl_plugin *plugin)
00199 {
00200
00201 cpl_recipe *recipe = (cpl_recipe *) plugin;
00202 if(recipe->parameters == NULL) {
00203 return 1;
00204 }
00205 if(recipe->frames == NULL) {
00206 return 1;
00207 }
00208
00209 return sinfo_utl_illumcorr(recipe->parameters, recipe->frames);
00210
00211 }
00212
00213 static int
00214 sinfo_utl_illumcorr_destroy(cpl_plugin *plugin)
00215 {
00216 cpl_recipe *recipe = (cpl_recipe *) plugin;
00217
00218
00219
00220
00221
00222
00223 cpl_parameterlist_delete(recipe->parameters);
00224
00225 return 0;
00226
00227 }
00228
00229
00230 int
00231 cpl_plugin_get_info(cpl_pluginlist *list)
00232 {
00233
00234 cpl_recipe *recipe = cpl_calloc(1, sizeof *recipe);
00235 cpl_plugin *plugin = &recipe->interface;
00236
00237 strcpy(sinfo_utl_illumcorr_description,sinfo_utl_illumcorr_description1);
00238 strcat(sinfo_utl_illumcorr_description,sinfo_utl_illumcorr_description2);
00239 strcat(sinfo_utl_illumcorr_description,sinfo_utl_illumcorr_description3);
00240 strcat(sinfo_utl_illumcorr_description,sinfo_utl_illumcorr_description4);
00241
00242 cpl_plugin_init(plugin,
00243 CPL_PLUGIN_API,
00244 SINFONI_BINARY_VERSION,
00245 CPL_PLUGIN_TYPE_RECIPE,
00246 "sinfo_utl_illumcorr",
00247 "Object data reduction",
00248 sinfo_utl_illumcorr_description,
00249 "Juha Reunanen",
00250 "reunanen@strw.leidenuniv.nl",
00251 sinfo_get_license(),
00252 sinfo_utl_illumcorr_create,
00253 sinfo_utl_illumcorr_exec,
00254 sinfo_utl_illumcorr_destroy);
00255
00256 cpl_pluginlist_append(list, plugin);
00257
00258 return 0;
00259
00260 }
00261
00262
00263
00264
00265
00266 static int
00267 sinfo_utl_illumcorr(cpl_parameterlist *config, cpl_frameset *set)
00268 {
00269 char outname[FILE_NAME_SZ];
00270
00271 int i=0;
00272 int k=0;
00273
00274 int ind=0;
00275 int nsky=0;
00276 int nobj=0;
00277 int ncdb=0;
00278 int nstk=0;
00279
00280 cpl_frameset * obj_set=NULL;
00281 cpl_frameset * sky_set=NULL;
00282 cpl_frameset * cdb_set=NULL;
00283 cpl_frameset * wrk_set=NULL;
00284 cpl_frameset * stk_set=NULL;
00285 cpl_frame * sky_frm=NULL;
00286
00287 cpl_frame * dup_frm=NULL;
00288 cpl_frame * cdb_frm=NULL;
00289 cpl_frame * wrk_frm=NULL;
00290 cpl_frameset * ref_set=NULL;
00291
00292 cpl_frame * dark_frm=NULL;
00293
00294 fake* fk;
00295
00296
00297 cpl_image * ima1=NULL ;
00298 cpl_image * ima2=NULL ;
00299 cpl_image * resima=NULL ;
00300 cpl_propertylist * plist=NULL ;
00301 cpl_frame * product_frame=NULL;
00302 const char *name_i=NULL;
00303
00304
00305 sinfo_msg("Welcome to SINFONI Pipeline release %d.%d.%d",
00306 SINFONI_MAJOR_VERSION,SINFONI_MINOR_VERSION,SINFONI_MICRO_VERSION);
00307
00308 if(sinfo_dfs_set_groups(set)) {
00309 sinfo_msg_error("Cannot identify RAW and CALIB frames") ;
00310 return -1;
00311 }
00312
00313 dark_frm = cpl_frameset_find(set,PRO_MASTER_DARK);
00314 if (dark_frm == NULL) {
00315 sinfo_msg_error("Cannot find dark frame") ;
00316 return (-1);
00317 }
00318
00319 ref_set=cpl_frameset_duplicate(set);
00320
00321 obj_set=cpl_frameset_new();
00322 sky_set=cpl_frameset_new();
00323 cdb_set=cpl_frameset_new();
00324 fk = sinfo_fake_new();
00325
00326 sinfo_extract_obj_frames(set,obj_set);
00327 sinfo_extract_sky_frames(set,sky_set);
00328 sinfo_extract_mst_frames(set,cdb_set);
00329
00330 nobj=cpl_frameset_get_size(obj_set);
00331 nsky=cpl_frameset_get_size(sky_set);
00332 ncdb=cpl_frameset_get_size(cdb_set);
00333
00334 if ((nobj==0) && (nsky==0)) {
00335 sinfo_msg_error("Empty input set");
00336 cpl_frameset_delete(obj_set);
00337 cpl_frameset_delete(sky_set);
00338 cpl_frameset_delete(cdb_set);
00339 cpl_frameset_delete(ref_set);
00340 sinfo_fake_delete(&fk);
00341 return (-1);
00342 }
00343
00344
00345
00346
00347
00348
00349 if ( nsky != 0) {
00350 if( (sky_frm = sinfo_get_dummy_object(sky_set)) == NULL) {
00351 sinfo_msg_error("Problem to get dummy frame");
00352 cpl_frameset_delete(obj_set);
00353 cpl_frameset_delete(sky_set);
00354 cpl_frameset_delete(cdb_set);
00355 cpl_frameset_delete(ref_set);
00356 sinfo_fake_delete(&fk);
00357 return (-1);
00358 }
00359 }
00360 else {
00361 if( (sky_frm = sinfo_get_dummy_object(obj_set)) == NULL) {
00362 sinfo_msg_error("Problem to get dummy frame");
00363 cpl_frameset_delete(obj_set);
00364 cpl_frameset_delete(sky_set);
00365 cpl_frameset_delete(cdb_set);
00366 cpl_frameset_delete(ref_set);
00367 sinfo_fake_delete(&fk);
00368 return (-1);
00369 }
00370 }
00371
00372
00373
00374
00375
00376 ima1 = cpl_image_load(cpl_frame_get_filename(sky_frm),CPL_TYPE_FLOAT,0,0);
00377 ima2 = cpl_image_load(cpl_frame_get_filename(dark_frm),CPL_TYPE_FLOAT,0,0);
00378 resima = cpl_image_subtract_create(ima1, ima2);
00379 plist=cpl_propertylist_load(cpl_frame_get_filename(sky_frm), 0);
00380 cpl_image_delete(ima1);
00381 cpl_image_delete(ima2);
00382
00383 product_frame = cpl_frame_new();
00384 cpl_frame_set_filename(product_frame, "out_fake_object2.fits") ;
00385 cpl_frame_set_tag(product_frame, "OBJECT_NODDING") ;
00386 cpl_frame_set_type(product_frame, CPL_FRAME_TYPE_IMAGE) ;
00387 cpl_frame_set_group(product_frame, CPL_FRAME_GROUP_RAW) ;
00388
00389 cpl_propertylist_erase_regexp(plist, "^ESO PRO CATG",0);
00390
00391 cpl_image_save(resima, "out_fake_object2.fits", CPL_BPP_DEFAULT, plist,
00392 CPL_IO_DEFAULT) ;
00393 cpl_propertylist_delete(plist) ;
00394 cpl_image_delete(resima) ;
00395
00396
00397
00398
00399
00400 wrk_set=cpl_frameset_new();
00401
00402 dup_frm=cpl_frame_duplicate(product_frame);
00403 cpl_frame_set_tag (dup_frm, "OBJECT_NODDING");
00404 cpl_frame_set_group (dup_frm ,CPL_FRAME_GROUP_RAW);
00405 cpl_frameset_insert(wrk_set,dup_frm);
00406
00407
00408 for(k=0;k<ncdb;k++) {
00409 cdb_frm=cpl_frameset_get_frame(cdb_set,k);
00410 dup_frm=cpl_frame_duplicate(cdb_frm);
00411 cpl_frameset_insert(wrk_set,dup_frm);
00412 }
00413
00414
00415
00416 sprintf(outname,"%s%d%s","out_stack",i,".fits");
00417 if(-1 == sinfo_new_stack_frames(config,wrk_set,
00418 PRO_OBJECT_NODDING_STACKED,i,fk,cpl_func)) {
00419
00420 cpl_frameset_delete(wrk_set);
00421
00422 cpl_frameset_delete(obj_set);
00423 cpl_frameset_delete(sky_set);
00424 cpl_frameset_delete(cdb_set);
00425 cpl_frameset_delete(ref_set);
00426 sinfo_fake_delete(&fk);
00427 return -1;
00428 }
00429
00430 stk_set=cpl_frameset_new();
00431 sinfo_contains_frames_kind(wrk_set,stk_set,PRO_STACKED);
00432 nstk=cpl_frameset_get_size(stk_set);
00433
00434 for(k=0;k<nstk;k++) {
00435 wrk_frm=cpl_frameset_get_frame(stk_set,k);
00436 dup_frm = cpl_frame_duplicate(wrk_frm);
00437 cpl_frameset_insert(set,dup_frm);
00438 }
00439 cpl_frameset_delete(stk_set);
00440 cpl_frameset_delete(wrk_set);
00441
00442 sinfo_msg("------------------------------") ;
00443 sinfo_msg("CREATING SKY CUBE");
00444 sinfo_msg("------------------------------") ;
00445
00446
00447 if ( -1 == (ind = sinfo_new_objnod(cpl_func,config, set, PRO_COADD_OBJ ) ) ) {
00448 sinfo_msg_error("NODDING SCIENCE FRAMES no. %d\n", ind) ;
00449 cpl_frameset_delete(obj_set);
00450 cpl_frameset_delete(sky_set);
00451 cpl_frameset_delete(cdb_set);
00452 cpl_frameset_delete(ref_set);
00453 sinfo_fake_delete(&fk);
00454
00455 return (-1);
00456 }
00457 sinfo_msg("------------------------------") ;
00458 sinfo_msg("CREATED SKY CUBE");
00459 sinfo_msg("------------------------------") ;
00460
00461
00462 stk_set=cpl_frameset_new();
00463 sinfo_contains_frames_kind(set, stk_set, PRO_OBS_OBJ);
00464 nstk=cpl_frameset_get_size(stk_set);
00465
00466 wrk_frm=cpl_frameset_get_frame(stk_set,0);
00467 name_i = cpl_frame_get_filename(wrk_frm);
00468
00469
00470
00471 cpl_frameset_delete(obj_set);
00472 cpl_frameset_delete(sky_set);
00473 cpl_frameset_delete(cdb_set);
00474 cpl_frameset_delete(ref_set);
00475 sinfo_fake_delete(&fk);
00476 cpl_frame_delete(sky_frm);
00477 create_illumcorr (cpl_func, config, set, name_i);
00478
00479 return (0);
00480
00481 }
00482
00483
00484 static cpl_frame*
00485 sinfo_get_dummy_object(cpl_frameset* obj_set)
00486 {
00487
00488 cpl_imagelist* obj_list=NULL;
00489 cpl_image* fake_object=NULL;
00490 char filename[FILE_NAME_SZ];
00491 cpl_frame* frame=NULL;
00492 cpl_frame* object_frame=NULL;
00493
00494 cpl_propertylist* plist=NULL;
00495
00496 obj_list = cpl_imagelist_load_frameset(obj_set,CPL_TYPE_FLOAT,0,0);
00497 fake_object = cpl_imagelist_collapse_median_create(obj_list);
00498
00499 frame = cpl_frameset_get_frame(obj_set,0);
00500 strcpy(filename,cpl_frame_get_filename(frame));
00501
00502 if ((cpl_error_code)((plist = cpl_propertylist_load(filename, 0)) == NULL)) {
00503 sinfo_msg_error("getting header from reference ima frame %s",filename);
00504 cpl_propertylist_delete(plist) ;
00505 return NULL ;
00506 }
00507
00508 if (cpl_propertylist_contains(plist, KEY_NAME_DPR_TYPE)) {
00509 cpl_propertylist_set_string(plist, KEY_NAME_DPR_TYPE, "OBJECT");
00510 } else {
00511 cpl_propertylist_append_string(plist, KEY_NAME_DPR_TYPE,"OBJECT") ;
00512 }
00513
00514 if (cpl_image_save(fake_object, "out_fake_object.fits", CPL_BPP_DEFAULT,
00515 plist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00516 sinfo_msg_error("Cannot save the product %s","out_fake_object.fits");
00517 cpl_propertylist_delete(plist) ;
00518 return NULL ;
00519 }
00520 cpl_propertylist_delete(plist);
00521
00522 object_frame = cpl_frame_new() ;
00523 cpl_frame_set_filename(object_frame, "out_fake_object.fits") ;
00524 cpl_frame_set_tag(object_frame, "OBJECT") ;
00525 cpl_frame_set_type(object_frame, CPL_FRAME_TYPE_IMAGE);
00526
00527
00528
00529 cpl_frame_set_level(object_frame, CPL_FRAME_LEVEL_FINAL);
00530 cpl_image_delete(fake_object);
00531 cpl_imagelist_delete(obj_list);
00532
00533 return object_frame;
00534 }
00535
00536 static void
00537 sinfo_illumcorr_config_add (cpl_parameterlist *list)
00538 {
00539 cpl_parameter *p;
00540
00541 if (!list) {
00542 return;
00543 }
00544
00545 p = cpl_parameter_new_range("sinfoni.illumcorr.spec_bin",
00546 CPL_TYPE_INT,
00547 "Number of spectral planes to be combined ",
00548 "sinfoni.illumcorr",
00549 100, 1, 200);
00550 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-spec_bin");
00551 cpl_parameterlist_append(list, p);
00552
00553 p = cpl_parameter_new_value("sinfoni.illumcorr.min_flux",
00554 CPL_TYPE_DOUBLE,
00555 "Minimum flux in each spectral bin ",
00556 "sinfoni.illumcorr",
00557 0.0);
00558 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-min_flux");
00559 cpl_parameterlist_append(list, p);
00560
00561 p = cpl_parameter_new_value("sinfoni.illumcorr.center_bins",
00562 CPL_TYPE_BOOL,
00563 "Center the spectral bins at prominent "
00564 "emission features ",
00565 "sinfoni.illumcorr",
00566 FALSE);
00567 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-center_bins");
00568 cpl_parameterlist_append(list, p);
00569
00570 p = cpl_parameter_new_enum("sinfoni.illumcorr.order",
00571 CPL_TYPE_INT,
00572 "The order of the polynomial to be fitted "
00573 "for each slitlet",
00574 "sinfoni.illumcorr",
00575 0,
00576 2,0,1);
00577 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-order");
00578 cpl_parameterlist_append(list, p);
00579
00580 p = cpl_parameter_new_value("sinfoni.illumcorr.sigma",
00581 CPL_TYPE_DOUBLE,
00582 "Reject n-sigma deviant pixels on each slitlet ",
00583 "sinfoni.illumcorr",
00584 3.0);
00585 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-sigma");
00586 cpl_parameterlist_append(list, p);
00587
00588 p = cpl_parameter_new_value("sinfoni.illumcorr.iterations",
00589 CPL_TYPE_INT,
00590 "Number of sigma rejection iterations to run ",
00591 "sinfoni.illumcorr",
00592 3);
00593 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-iter");
00594 cpl_parameterlist_append(list, p);
00595
00596 p = cpl_parameter_new_range("sinfoni.illumcorr.llx",
00597 CPL_TYPE_INT,
00598 "Reference region coordinates ",
00599 "sinfoni.illumcorr",
00600 8, 0, 63);
00601 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-llx");
00602 cpl_parameterlist_append(list, p);
00603
00604 p = cpl_parameter_new_range("sinfoni.illumcorr.lly",
00605 CPL_TYPE_INT,
00606 "Reference region coordinates ",
00607 "sinfoni.illumcorr",
00608 33, 0, 63);
00609 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-lly");
00610 cpl_parameterlist_append(list, p);
00611
00612 p = cpl_parameter_new_range("sinfoni.illumcorr.urx",
00613 CPL_TYPE_INT,
00614 "Reference region coordinates ",
00615 "sinfoni.illumcorr",
00616 60, 0, 63);
00617 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-urx");
00618 cpl_parameterlist_append(list, p);
00619
00620 p = cpl_parameter_new_range("sinfoni.illumcorr.ury",
00621 CPL_TYPE_INT,
00622 "Reference region coordinates ",
00623 "sinfoni.illumcorr",
00624 36, 0, 63);
00625 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-ury");
00626 cpl_parameterlist_append(list, p);
00627
00628 p = cpl_parameter_new_enum("sinfoni.illumcorr.smooth0",
00629 CPL_TYPE_INT,
00630 "Smooth zeroth order terms by fitting "
00631 "with polynomial (1),"
00632 "with median filter (2) or not (0) ",
00633 "sinfoni.illumcorr",
00634 0,
00635 3, 0, 1, 2);
00636 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-smooth0");
00637 cpl_parameterlist_append(list, p);
00638
00639 p = cpl_parameter_new_value("sinfoni.illumcorr.smooth0_order",
00640 CPL_TYPE_INT,
00641 "Order of the smoothing polynomial for 0th term",
00642 "sinfoni.illumcorr",
00643 2);
00644 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-smooth_order0");
00645 cpl_parameterlist_append(list, p);
00646
00647 p = cpl_parameter_new_range("sinfoni.illumcorr.smooth0_size",
00648 CPL_TYPE_INT,
00649 "Size of the median filter for 0th "
00650 "order smoothing ",
00651 "sinfoni.illumcorr",
00652 51,3, 301);
00653 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-smooth0_size");
00654 cpl_parameterlist_append(list, p);
00655
00656 p = cpl_parameter_new_value("sinfoni.illumcorr.smooth1",
00657 CPL_TYPE_BOOL,
00658 "Smooth higher (>0) order polynomials ",
00659 "sinfoni.illumcorr",
00660 TRUE);
00661 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-smooth");
00662 cpl_parameterlist_append(list, p);
00663
00664 p = cpl_parameter_new_value("sinfoni.illumcorr.smooth1_order",
00665 CPL_TYPE_INT,
00666 "Smoothing order for higher terms ",
00667 "sinfoni.illumcorr",
00668 2);
00669 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-smooth_order");
00670 cpl_parameterlist_append(list, p);
00671
00672 p = cpl_parameter_new_value("sinfoni.illumcorr.illumcorr_sigma",
00673 CPL_TYPE_DOUBLE,
00674 "Reject all fits for which the rms is "
00675 "illumcorr-sigma times bigger than the "
00676 "median rms in each spectral bin" ,
00677 "sinfoni.illumcorr",
00678 5.0);
00679 cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,
00680 "illumcorr-illumcorr_sigma");
00681 cpl_parameterlist_append(list, p);
00682
00683 }
00684
00685 static int
00686 create_illumcorr (const char* plugin_id,
00687 cpl_parameterlist *cpl_cfg,
00688 cpl_frameset* sof,
00689 const char *name_i)
00690 {
00691 cpl_parameter *p=NULL;
00692 double min_flux=0;
00693 double sigma=0;
00694 double il_sigma=0;
00695 int spec_bin=0;
00696 int _order=0;
00697 cpl_imagelist *sky=NULL;
00698 cpl_imagelist *binnedsky=NULL;
00699 cpl_imagelist *result=NULL;
00700 cpl_image *temp_image=NULL;
00701 cpl_image *temp_image2=NULL;
00702 int nplanes=0;
00703 int i=0;
00704 int j=0;
00705 int k=0;
00706 int kk=0;
00707 int n=0;
00708 int slitlet=0;
00709 int bin=0;
00710 double *median=NULL;
00711 double *pos=NULL;
00712 double temp=0;
00713 double temp2=0;
00714 double *inter_pos=NULL;
00715 double *inter_val=NULL;
00716 double *plane_pos=NULL;
00717 double *plane_val=NULL;
00718 int llx=0;
00719 int lly=0;
00720 int urx=0;
00721 int ury=0;
00722 int smooth_order=0;
00723 int iter=0;
00724 cpl_vector *row=NULL;
00725 cpl_vector *model=NULL;
00726 cpl_vector *xpos=NULL;
00727 cpl_vector *tempvector=NULL;
00728 cpl_vector *tempvector2=NULL;
00729 double mse=0.0;
00730 double stddev=0.0;
00731 cpl_polynomial*poly=NULL;
00732 cpl_polynomial *poly2=NULL;
00733 double *temparray=NULL;
00734 double *tempxarray=NULL;
00735 double * tempsarray=NULL;
00736 cpl_polynomial**coeffs=NULL;
00737 float *data=NULL;
00738 double *rms_values=NULL;
00739 double rms_array[32];
00740 int smooth=0;
00741 int smooth0=0;
00742 int smooth_order0=0;
00743 int smooth_size0=0;
00744 int center_bins = 0;
00745
00746 int *bin_start=NULL;
00747 int *bin_end=NULL;
00748 int z1=0;
00749 int z2=0;
00750 int nbins=0;
00751
00752 FILE *dumpfile=NULL;
00753
00754 int order[32];
00755 int ok[64];
00756 int nbad=0;
00757
00758
00759
00760
00761
00762 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.spec_bin");
00763 spec_bin = cpl_parameter_get_int(p);
00764 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.min_flux");
00765 min_flux = cpl_parameter_get_double(p);
00766 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.order");
00767 _order = cpl_parameter_get_int(p);
00768 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.sigma");
00769 sigma = cpl_parameter_get_double(p);
00770 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.llx");
00771 llx = cpl_parameter_get_int(p);
00772 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.lly");
00773 lly = cpl_parameter_get_int(p);
00774 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.urx");
00775 urx = cpl_parameter_get_int(p);
00776 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.ury");
00777 ury = cpl_parameter_get_int(p);
00778 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.illumcorr_sigma");
00779 il_sigma = cpl_parameter_get_double(p);
00780
00781 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.smooth0");
00782 smooth0 = cpl_parameter_get_int (p);
00783 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.smooth0_order");
00784 smooth_order0 = cpl_parameter_get_int (p);
00785 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.smooth0_size");
00786 smooth_size0 = cpl_parameter_get_int (p);
00787
00788 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.smooth1");
00789 smooth = cpl_parameter_get_bool (p);
00790 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.smooth1_order");
00791 smooth_order = cpl_parameter_get_int (p);
00792
00793 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.iterations");
00794 iter = cpl_parameter_get_int (p);
00795
00796 p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.center_bins");
00797 center_bins = cpl_parameter_get_bool (p);
00798
00799
00800
00801
00802
00803
00804 sky = cpl_imagelist_load(name_i, CPL_TYPE_FLOAT, 0);
00805 nplanes = cpl_imagelist_get_size(sky);
00806
00807
00808
00809 z1 = 0;
00810 z2=nplanes -1;
00811 while (z1<nplanes
00812 && isnan(cpl_image_get_mean_window(cpl_imagelist_get(sky, z1),
00813 llx, lly, urx, ury)))
00814 z1++;
00815 while (z2>=0
00816 && isnan(cpl_image_get_mean_window(cpl_imagelist_get(sky, z2),
00817 llx, lly, urx, ury)))
00818 z2--;
00819 z1 += 2;
00820 z2 -= 2;
00821 if (z1>=nplanes || z2 <0 || z2<=z1) {
00822 sinfo_msg_error ("Start z = %d, end z = %d", z1, z2);
00823 cpl_imagelist_delete (sky);
00824 return (-1);
00825 }
00826 sinfo_msg ("Start z = %d, end z = %d", z1, z2);
00827
00828 binnedsky = cpl_imagelist_new ();
00829 median = (double*) cpl_calloc(nplanes/spec_bin, sizeof(double));
00830 pos = (double*) cpl_calloc(nplanes/spec_bin, sizeof(double));
00831 temparray = (double*) cpl_calloc(64, sizeof(double));
00832 tempxarray= (double*) cpl_calloc(64, sizeof(double));
00833 tempsarray= (double*) cpl_calloc (nplanes/spec_bin, sizeof(double));
00834 plane_pos = (double*) cpl_calloc (nplanes/spec_bin, sizeof(double));
00835 plane_val = (double*) cpl_calloc (nplanes/spec_bin, sizeof(double));
00836 coeffs = (cpl_polynomial**) cpl_calloc(32*(nplanes/spec_bin),
00837 sizeof(cpl_polynomial*));
00838 rms_values= (double*) cpl_calloc (32*(nplanes/spec_bin), sizeof (double));
00839 inter_pos = (double*) cpl_calloc (nplanes, sizeof(double));
00840 inter_val = (double*) cpl_calloc (nplanes, sizeof(double));
00841
00842 model = cpl_vector_new(64);
00843 xpos = cpl_vector_new(64);
00844
00845 for (i=0; i<64; i++)
00846 cpl_vector_set(xpos, i, (double)(i)-((double)urx-(double)llx)/2.0);
00847 for (i=0; i<nplanes; i++)
00848 inter_pos[i] = (double)i;
00849
00850
00851
00852
00853
00854
00855 for (i=0; i<32; i++)
00856 order[i] = _order;
00857
00858
00859 if (center_bins == 1) {
00860 sinfo_msg("Using centering on emission features\n");
00861 nbins = sinfo_illumcorr_create_bins (sky,llx, lly, urx, ury,
00862 spec_bin, min_flux,
00863 &bin_start, &bin_end,
00864 z1, z2);
00865 }
00866 else {
00867 sinfo_msg("Using simple spectral binning - "
00868 "not centering on emission features\n");
00869 nbins = (z2-z1)/spec_bin;
00870 bin_start = (int*)cpl_calloc(nbins+1, sizeof(int));
00871 bin_end = (int*)cpl_calloc(nbins+1, sizeof(int));
00872 for (i = 0; i<nbins; i++) {
00873 bin_start[i] = z1+i*spec_bin;
00874 bin_end[i] = z1+(i+1)*spec_bin - 1;
00875 }
00876 if (bin_end[nbins-1]<z2-spec_bin/10) {
00877 bin_start[nbins] = bin_end[nbins-1]+1;
00878 bin_end[nbins] = z2;
00879 nbins++;
00880 }
00881 }
00882
00883
00884
00885
00886
00887
00888
00889 sinfo_msg("Binning the cube and calculating statistics\n");
00890 for (i=0; i<nbins; i++) {
00891 temp_image = cpl_image_duplicate(cpl_imagelist_get(sky, bin_start[i]));
00892 median[i] = sinfo_image_get_median_window (temp_image, llx, lly, urx, ury);
00893 pos[i] = median[i] * (double)bin_start[i];
00894 cpl_imagelist_set (binnedsky, temp_image, i);
00895 for (j=bin_start[i]+1; j<bin_end[i]; j++) {
00896 temp_image2 = cpl_imagelist_get (sky, j);
00897 cpl_image_add (temp_image, temp_image2);
00898 temp = sinfo_image_get_median_window (temp_image2, llx, lly, urx, ury);
00899 median[i] = median[i] + temp;
00900 pos[i] = pos[i] + temp*(double)j;
00901 }
00902 temp2 =(double)(bin_end[i]-bin_start[i]+1);
00903 cpl_image_divide_scalar (temp_image, temp2);
00904 pos[i] = pos[i]/median[i];
00905 median[i] = median[i] / temp2;
00906 sinfo_msg_debug("median image=%g at %g",median[i], pos[i]);
00907 }
00908
00909 sinfo_msg("Fitting slitlets between x=%d - x=%d\n", llx, urx);
00910 sinfo_msg("Fitting order %d\n", _order);
00911 for (k=0;k<nbins; k++) {
00912 if (median[k]>min_flux) {
00913 for (j=0; j<32; j++) {
00914 row=cpl_vector_new_from_image_row(cpl_imagelist_get(binnedsky,k),2*j+1);
00915 n = 0;
00916 for (i=llx; i<=urx; i++) {
00917 if (!isnan(cpl_vector_get(row, i))) {
00918 ok[i] = 1;
00919 temparray[n] = cpl_vector_get(row, i);
00920 tempxarray[n] = cpl_vector_get(xpos, i);
00921 n++;
00922 }
00923 else
00924 ok[i] = 0;
00925 }
00926
00927
00928 if (n>20) {
00929 tempvector = cpl_vector_wrap (n, temparray);
00930 tempvector2= cpl_vector_wrap (n, tempxarray);
00931 poly = cpl_polynomial_fit_1d_create (tempvector2, tempvector,
00932 order[j], &mse);
00933
00934 if (poly == NULL)
00935 sinfo_msg("Fitting failed (plane %d, row %d): %s",
00936 k, j, (char* ) cpl_error_get_message());
00937 else {
00938 if (sigma>0 && iter>0) {
00939 for (kk = 0; kk<iter; kk++) {
00940 cpl_vector_fill_polynomial (model, poly, 0.0, 1.0);
00941 cpl_vector_subtract (model, row);
00942
00943
00944 n = 0;
00945 for (i=llx; i<=urx; i++)
00946 if (ok[i] == 1)
00947 temparray[n++] = cpl_vector_get(model, i);
00948 stddev = cpl_vector_get_stdev(tempvector);
00949
00950 for (i=llx; i<=urx; i++)
00951 if (ok[i] == 1)
00952 if (fabs(cpl_vector_get(model, i))>(stddev*sigma))
00953 ok[i] = 0;
00954
00955
00956 n = 0;
00957 for (i=llx; i<=urx; i++) {
00958 if (ok[i] == 1) {
00959 temparray[n] = cpl_vector_get(row, i);
00960 tempxarray[n] = cpl_vector_get(xpos, i);
00961 n++;
00962 }
00963 }
00964 cpl_polynomial_delete(poly);
00965 if (n>20) {
00966 cpl_vector_unwrap (tempvector);
00967 cpl_vector_unwrap (tempvector2);
00968 tempvector = cpl_vector_wrap (n, temparray);
00969 tempvector2= cpl_vector_wrap (n, tempxarray);
00970 stddev = cpl_vector_get_stdev(tempvector);
00971
00972 poly = cpl_polynomial_fit_1d_create (tempvector2,
00973 tempvector,
00974 order[j], &mse);
00975 if (poly == NULL)
00976 break;
00977 }
00978 else {
00979 poly = NULL;
00980 break;
00981 }
00982
00983 }
00984 }
00985
00986 if (poly!=NULL) {
00987 coeffs[j*nbins+k] = poly;
00988 rms_values[j*nbins+k] = sqrt(mse);
00989 }
00990 else
00991 coeffs[j*nbins+k] = NULL;
00992 }
00993 cpl_vector_unwrap (tempvector);
00994 cpl_vector_unwrap (tempvector2);
00995 }
00996 cpl_vector_delete(row);
00997 }
00998 }
00999 }
01000
01001
01002
01003
01004 sinfo_msg("Writing dat out_illum.dat\n");
01005 dumpfile = fopen ("out_illum.dat","w");
01006 fprintf (dumpfile, "# slitlet, pos, median, rms, coeff0, coeff1...\n");
01007 for (slitlet = 0; slitlet<32; slitlet++)
01008 for (bin=0; bin<nbins; bin++) {
01009 poly = coeffs[slitlet*nbins+bin];
01010 if (poly != NULL) {
01011 fprintf (dumpfile, "%d %f %f %f ",slitlet, pos[bin],
01012 median[bin],
01013 rms_values[slitlet*nbins+bin]);
01014 for (i=0; i<(cpl_polynomial_get_degree(poly)+1); i++) {
01015 temp = cpl_polynomial_get_coeff (poly, &i);
01016 fprintf (dumpfile, "%f ", temp);
01017 }
01018 fprintf (dumpfile, "\n");
01019 }
01020 }
01021 fclose (dumpfile);
01022
01023
01024
01025
01026
01027
01028 sinfo_msg("Removing poor fits - factor %f", il_sigma);
01029 n = 0;
01030 i = 0;
01031 nbad=0;
01032 sinfo_msg("max loop over bin =%d",nbins);
01033 for (bin=0; bin<nbins; bin++) {
01034 k = 0;
01035 for (slitlet=0; slitlet<32; slitlet++)
01036 if (coeffs[slitlet*nbins+bin] != NULL)
01037 rms_array[k++] = rms_values[slitlet*nbins+bin];
01038 if (k>0) {
01039
01040
01041
01042 tempvector = cpl_vector_wrap (k, &rms_array[0]);
01043 temp = cpl_vector_get_median (tempvector);
01044 sinfo_msg("median temp=%g",temp);
01045 cpl_vector_unwrap (tempvector);
01046 for (slitlet=0; slitlet<32; slitlet++)
01047 if (coeffs[slitlet*nbins+bin] != NULL) {
01048 i++;
01049 if (rms_values[slitlet*nbins+bin]>(il_sigma*temp)) {
01050 cpl_polynomial_delete(coeffs[slitlet*nbins+bin]);
01051 coeffs[slitlet*nbins+bin] = NULL;
01052 n++;
01053 }
01054 } else {
01055 nbad++;
01056 }
01057
01058 }
01059 }
01060 sinfo_msg("Removed %d poor fits out of %d. Bad coeffs=%d", n,i,nbad);
01061
01062 if(smooth0 == 1) {
01063 sinfo_msg("Smoothing zeroth terms (order %d)", smooth_order0);
01064
01065
01066
01067 for (slitlet = 0; slitlet<32; slitlet++) {
01068 k = 0;
01069 for (bin=0; bin<nbins; bin++) {
01070 if (coeffs[slitlet*nbins+bin] != NULL) {
01071 poly = coeffs[slitlet*nbins+bin];
01072 i = 0;
01073 temp = cpl_polynomial_get_coeff (poly, &i);
01074 plane_pos[k] = pos[bin];
01075 plane_val[k] = temp/median[bin];
01076 k++;
01077 }
01078 }
01079 if (k>0) {
01080 tempvector = cpl_vector_wrap (k, plane_pos);
01081 tempvector2= cpl_vector_wrap (k, plane_val);
01082 poly2 = cpl_polynomial_fit_1d_create (tempvector, tempvector2,
01083 smooth_order0, &mse);
01084 cpl_vector_unwrap (tempvector);
01085 cpl_vector_unwrap (tempvector2);
01086 for (bin=0; bin<nbins; bin++) {
01087 if (coeffs[slitlet*nbins+bin] != NULL) {
01088 poly = coeffs[slitlet*nbins+bin];
01089 i = 0;
01090 temp2 = cpl_polynomial_eval_1d (poly2, pos[bin], NULL);
01091 cpl_polynomial_set_coeff (poly, &i, temp2*median[bin]);
01092 }
01093 }
01094 cpl_polynomial_delete(poly2);
01095 }
01096 else
01097 sinfo_msg_warning ("Not enough data points in slitlet %d", slitlet);
01098 }
01099 }
01100 else if (smooth0 == 2) {
01101 sinfo_msg("Smoothing zeroth terms (median filter size %d)", smooth_size0);
01102 smooth_size0 = smooth_size0/2;
01103 if (smooth_size0 <= 0) smooth_size0 = 1;
01104 for (slitlet = 0; slitlet<32; slitlet++) {
01105 k = 0;
01106 for (bin=0; bin<nbins; bin++) {
01107 if (coeffs[slitlet*nbins+bin] != NULL) {
01108 poly = coeffs[slitlet*nbins+bin];
01109 i = 0;
01110 temp = cpl_polynomial_get_coeff (poly, &i);
01111
01112 plane_val[k] = temp/median[bin];
01113 k++;
01114 }
01115 }
01116 if (k>0) {
01117 tempvector = cpl_vector_wrap (k, plane_val);
01118 tempvector2= sinfo_juha_vector_filter_median_create (tempvector,
01119 smooth_size0);
01120 cpl_vector_unwrap (tempvector);
01121 kk = 0;
01122 for (bin=0; bin<nbins; bin++) {
01123 if (coeffs[slitlet*nbins+bin] != NULL) {
01124 poly = coeffs[slitlet*nbins+bin];
01125 i = 0;
01126 cpl_polynomial_set_coeff(poly, &i, cpl_vector_get(tempvector2, kk++)
01127 *median[bin]);
01128 }
01129 }
01130 cpl_vector_delete (tempvector2);
01131 }
01132 }
01133 }
01134
01135 if(smooth == 1) {
01136 sinfo_msg("Smoothing higher terms (with order %d)", smooth_order);
01137 for (slitlet = 0; slitlet<32; slitlet++) {
01138 if (order[slitlet]>0) {
01139 for (j=1; j<=order[slitlet]; j++) {
01140 k = 0;
01141 for (bin=0; bin<nbins; bin++) {
01142 if (coeffs[slitlet*nbins+bin] != NULL) {
01143 poly = coeffs[slitlet*nbins+bin];
01144 i = 0;
01145
01146 temp2 = cpl_polynomial_get_coeff (poly, &j);
01147 plane_pos[k] = pos[bin];
01148 plane_val[k] = temp2/median[bin];
01149 k++;
01150 }
01151 }
01152 if (k>0) {
01153 tempvector = cpl_vector_wrap (k, plane_pos);
01154 tempvector2= cpl_vector_wrap (k, plane_val);
01155 poly2 = cpl_polynomial_fit_1d_create (tempvector, tempvector2,
01156 smooth_order, &mse);
01157 cpl_vector_unwrap (tempvector);
01158 cpl_vector_unwrap (tempvector2);
01159 for (bin=0; bin<nbins; bin++) {
01160 if (coeffs[slitlet*nbins+bin] != NULL) {
01161 poly = coeffs[slitlet*nbins+bin];
01162 i = 0;
01163
01164 temp2 = cpl_polynomial_eval_1d (poly2, pos[bin], NULL);
01165 cpl_polynomial_set_coeff (poly, &j, temp2*median[bin]);
01166 }
01167 }
01168 cpl_polynomial_delete(poly2);
01169 }
01170 else
01171 sinfo_msg_warning ("Not enough data points in slitlet %d\n",
01172 slitlet);
01173 }
01174 }
01175 }
01176 }
01177
01178 sinfo_msg("Creating cube for illumination correction\n");
01179 result = cpl_imagelist_new ();
01180 for (i=0; i<nplanes; i++) {
01181 temp_image = cpl_image_new (64, 64, CPL_TYPE_FLOAT);
01182 cpl_imagelist_set (result, temp_image, i);
01183 }
01184
01185 sinfo_msg("ok1 nplanes=%d spec_bin=%d",nplanes,spec_bin);
01186 if ( nbins>5) {
01187 sinfo_msg("Interpolating\n");
01188 for (slitlet = 0; slitlet<32; slitlet++) {
01189 for (i=0; i<64; i++) {
01190 k = 0;
01191 for (bin=0; bin<nbins; bin++) {
01192 if (coeffs[slitlet*nbins+bin] != NULL) {
01193 plane_pos[k] = pos[bin];
01194 plane_val[k] = cpl_polynomial_eval_1d(coeffs[slitlet*nbins+bin],
01195 cpl_vector_get(xpos, i),NULL)/
01196 median[bin];
01197 k++;
01198 }
01199 }
01200
01201 if (k>3) {
01202 sinfo_juha_function1d_natural_spline (plane_pos, plane_val, k,
01203 &inter_pos[(int)plane_pos[0]],
01204 &inter_val[(int)plane_pos[0]],
01205 (int)(plane_pos[k-1]-plane_pos[0]));
01206 for (j=0; j<=(int)plane_pos[0]; j++)
01207 inter_val[j] = inter_val[(int)plane_pos[0]+1];
01208 for (j=(int)plane_pos[k-1]-1; j<nplanes; j++)
01209 inter_val[j] = inter_val[(int)plane_pos[k-1]-2];
01210 for (k=0; k<nplanes; k++) {
01211 data = cpl_image_get_data_float(cpl_imagelist_get(result, k));
01212 data[i + (2*slitlet)*64] = inter_val[k];
01213 data[i + (2*slitlet+1)*64] = inter_val[k];
01214
01215 }
01216 }
01217 else
01218 sinfo_msg_warning ("Too few points %d\n", slitlet);
01219 }
01220 }
01221 }
01222 else if (nbins==1) {
01223 sinfo_msg("Filling the illumination cube\n");
01224 for (slitlet = 0; slitlet<32; slitlet++) {
01225 for (i=0; i<64; i++) {
01226 if (coeffs[slitlet] != NULL) {
01227 temp = cpl_polynomial_eval_1d(coeffs[slitlet],
01228 cpl_vector_get(xpos, i),NULL)/median[0];
01229 for (k=0; k<nplanes; k++) {
01230 data = cpl_image_get_data_float(cpl_imagelist_get(result, k));
01231 data[i + (2*slitlet)*64] = temp;
01232 data[i + (2*slitlet+1)*64] = temp;
01233 sinfo_msg("temp=%g",temp);
01234 }
01235 }
01236 }
01237 }
01238 } else {
01239
01240 }
01241 sinfo_msg("ok2");
01242
01243
01244 sinfo_msg("Writing ima out_illum.fits\n");
01245
01246
01247
01248
01249
01250
01251
01252
01253 sinfo_pro_save_ims(result,sof,sof,"out_illum.fits",
01254 PRO_ILL_COR,NULL,plugin_id, cpl_cfg);
01255
01256
01257
01258
01259 sinfo_msg("Writing dat out_illum2.dat\n");
01260 dumpfile = fopen ("out_illum2.dat","w");
01261 fprintf (dumpfile, "# slitlet, pos, median, rms, coeff0, coeff1...\n");
01262 for (slitlet = 0; slitlet<32; slitlet++)
01263 for (bin=0; bin<nbins; bin++) {
01264 poly = coeffs[slitlet*nbins+bin];
01265 if (poly != NULL) {
01266 fprintf (dumpfile, "%d %f %f %f ",slitlet, pos[bin],
01267 median[bin],
01268 rms_values[slitlet*nbins+bin]);
01269 for (i=0; i<(cpl_polynomial_get_degree(poly)+1); i++) {
01270 temp = cpl_polynomial_get_coeff (poly, &i);
01271 fprintf (dumpfile, "%f ", temp);
01272 }
01273 fprintf (dumpfile, "\n");
01274 }
01275 }
01276 fclose (dumpfile);
01277
01278
01279
01280
01281 for (i = 0; i<32*nbins; i++)
01282 if (coeffs[i] != NULL)
01283 cpl_polynomial_delete(coeffs[i]);
01284 cpl_imagelist_delete (sky);
01285 cpl_imagelist_delete (binnedsky);
01286 cpl_imagelist_delete (result);
01287 cpl_free (pos);
01288 cpl_free (median);
01289 cpl_free (temparray);
01290 cpl_free (tempxarray);
01291 cpl_free (tempsarray);
01292 cpl_free (coeffs);
01293 cpl_free (inter_pos);
01294 cpl_free (inter_val);
01295 cpl_free (plane_pos);
01296 cpl_free (plane_val);
01297 cpl_free (rms_values);
01298 cpl_vector_delete (xpos);
01299 cpl_vector_delete (model);
01300
01301 cpl_free (bin_start);
01302 cpl_free (bin_end);
01303
01304 return (1);
01305 }
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321 static int
01322 sinfo_illumcorr_create_bins (cpl_imagelist *sky,
01323 int llx, int lly, int urx, int ury,
01324 int spec_bin,
01325 double min_flux,
01326 int ** start,
01327 int ** end,
01328 int z1, int z2)
01329 {
01330 int temp_i=0;
01331 double testarray3[15];
01332 double temp_double=0;
01333 int i=0, j=0, k=0,kk=0,nplanes=0;
01334
01335 int norig = 0, nmerged = 0, ncont = 0, nline=0;
01336
01337 int *pos=NULL;
01338 int *x1=NULL;
01339 int *x2=NULL;
01340 int *x1b=NULL;
01341 int *x2b=NULL;
01342 int *s1=NULL;
01343 int *s2=NULL;
01344 double *flux=NULL;
01345 double *spec=NULL;
01346 double *spec_cont=NULL;
01347 double *spec_line=NULL;
01348
01349
01350 cpl_image *temp_image=NULL;
01351
01352 nplanes = cpl_imagelist_get_size(sky);
01353
01354 spec = (double*) cpl_calloc(nplanes, sizeof(double));
01355 spec_line = (double*) cpl_calloc(nplanes, sizeof(double));
01356 spec_cont = (double*) cpl_calloc(nplanes, sizeof(double));
01357
01358
01359 pos = (int*) cpl_calloc(nplanes, sizeof(int));
01360 flux = (double*) cpl_calloc(nplanes, sizeof(double));
01361 x1 = (int*) cpl_calloc(nplanes, sizeof(int));
01362 x2 = (int*) cpl_calloc(nplanes, sizeof(int));
01363 x1b = (int*) cpl_calloc(nplanes, sizeof(int));
01364 x2b = (int*) cpl_calloc(nplanes, sizeof(int));
01365
01366 for (i=z1; i<=z2; i++) {
01367 temp_image = cpl_imagelist_get(sky, i);
01368 spec[i] = sinfo_image_get_median_window (temp_image, llx, lly, urx, ury);
01369 }
01370 for (i=z1+7; i<=z2-7; i++) {
01371 k = 0;
01372 for (j=-7; j<=7; j++)
01373 if (!isnan(spec[i+j]))
01374 testarray3[k++] = spec[i+j];
01375 if (k>0) {
01376 sinfo_tools_sort_double (&testarray3[0], k);
01377 spec_cont[i] = testarray3[1];
01378 }
01379 else
01380 spec_cont[i] = 0./0.;
01381 }
01382
01383 sinfo_msg_debug("Calculating pure line flux at pos: "
01384 "original, continuum, line");
01385 for (i=z1; i<=z2; i++) {
01386 spec_line[i] = spec[i] - spec_cont[i];
01387 sinfo_msg_debug("Flux at %i = %g %g %g",
01388 i,spec[i],spec_cont[i], spec_line[i]);
01389 }
01390
01391
01392
01393
01394
01395 sinfo_msg ("Searching for peaks");
01396 temp_double = -10000.0;
01397 i = z1+2;
01398 while (i<=z2-2) {
01399 if (!isnan (spec_line[i])) {
01400 if (temp_double<spec_line[i]) {
01401 temp_i = i;
01402 temp_double = spec_line[i];
01403 }
01404 else {
01405
01406 if (temp_i == i-1 && spec_line[temp_i]>min_flux) {
01407 k = 0;
01408 for (j=-spec_bin/2; j<=spec_bin/2; j++)
01409 if (j+i>=0 && i+j<nplanes && isnan(spec[i+j])) {
01410 k = 1;
01411 break;
01412 }
01413 if (k==0) {
01414 pos[norig] = temp_i;
01415 flux[norig] = temp_double;
01416 x1[norig] = temp_i;
01417 x2[norig] = temp_i;
01418 temp_double = -10000.0;
01419 while (spec_line[i]<spec_line[i-1])
01420 i++;
01421 i--;
01422 norig++;
01423 }
01424 }
01425 }
01426 }
01427 i++;
01428 }
01429
01430
01431
01432
01433 sinfo_msg ("Merging emission features too close to each other");
01434 i = 0;
01435 while (i<norig) {
01436 if (flux[i] > 0.0) {
01437 j = i+1;
01438 while (j<norig
01439 && (x1[j]-x2[i]) <=spec_bin
01440 && flux[j]>0.0) {
01441 flux[j] = -100.0;
01442 x2[i] = x1[j];
01443 j++;
01444 nmerged++;
01445 }
01446 }
01447 i++;
01448 }
01449
01450 nline = norig - nmerged;
01451
01452 j = 0;
01453 for (i=0; i<norig; i++)
01454 if (flux[i]>0.0) {
01455 x1b[j] = x1[i] - spec_bin/2;
01456 x2b[j] = x2[i] + spec_bin/2;
01457 j++;
01458
01459 }
01460
01461 x1b[j] = nplanes+1;
01462
01463
01464
01465
01466
01467 j=0;
01468 i=z1;
01469 while (i<=z2) {
01470 if (!isnan (spec[i])) {
01471 if (x1b[j]-i < spec_bin) {
01472 i = x2b[j]+1;
01473 j++;
01474 }
01475 else {
01476 kk = 0;
01477 for (k=0; k<spec_bin; k++)
01478 if (spec[i+k]>min_flux)
01479 kk++;
01480 if (kk==spec_bin) {
01481 x1[ncont] = i;
01482 x2[ncont] = i+spec_bin-1;
01483 ncont++;
01484 i = i+spec_bin;
01485 }
01486 else
01487 i++;
01488 }
01489 }
01490 else
01491 i++;
01492 }
01493
01494 sinfo_msg ("Number of bins centered on emission features:");
01495 sinfo_msg (" %i - %i (merged) = %i", norig, nmerged, nline);
01496 sinfo_msg (" %i continuum bins", ncont);
01497
01498 s1 = (int*)cpl_calloc(norig-nmerged+ncont, sizeof(int));
01499 s2 = (int*)cpl_calloc(norig-nmerged+ncont, sizeof(int));
01500
01501
01502
01503
01504
01505 i=0;
01506 j=0;
01507 k=0;
01508 while (k<norig-nmerged+ncont) {
01509 if (i<norig && j<ncont && x1b[i]<x1[j]) {
01510 s1[k] = x1b[i];
01511 s2[k] = x2b[i];
01512 k++;
01513 i++;
01514 }
01515 else if (i<norig && j<ncont && x1b[i]>x1[j]) {
01516 s1[k] = x1[j];
01517 s2[k] = x2[j];
01518 k++;
01519 j++;
01520 }
01521 else if (i == norig) {
01522 s1[k] = x1[j];
01523 s2[k] = x2[j];
01524 k++;
01525 j++;
01526 }
01527 else if (j == ncont) {
01528 s1[k] = x1b[i];
01529 s2[k] = x2b[i];
01530 k++;
01531 i++;
01532 }
01533 else {
01534
01535 sinfo_msg_warning("Something went wrong when combining "
01536 "the bins %i and %i", i,j);
01537 break;
01538 }
01539 }
01540
01541 for (i=0; i<nline+ncont; i++)
01542 sinfo_msg_debug ("Bin start: %i end %i", s1[i], s2[i]);
01543
01544 *start = s1;
01545 *end = s2;
01546
01547 cpl_free (pos);
01548 cpl_free (x1);
01549 cpl_free (x2);
01550 cpl_free (x1b);
01551 cpl_free (x2b);
01552 cpl_free (flux);
01553 cpl_free (spec);
01554 cpl_free (spec_line);
01555 cpl_free (spec_cont);
01556
01557 return (nline+ncont);
01558 }
01559
01560
01561
01588
01589
01590 static int
01591 sinfo_juha_function1d_natural_spline(
01592 double * x,
01593 double * y,
01594 int len,
01595 double * splX,
01596 double * splY,
01597 int splLen
01598 )
01599 {
01600 int end;
01601 int loc, found;
01602 register int i, j, n;
01603 double * h;
01604 double * alpha;
01605 double * l,
01606 * mu,
01607 * z,
01608 * a,
01609 * b,
01610 * c,
01611 * d,
01612 v;
01613
01614 end = len - 1;
01615
01616 a = cpl_malloc(sizeof(double) * splLen * 9) ;
01617 b = a + len;
01618 c = b + len;
01619 d = c + len;
01620 h = d + len;
01621 l = h + len;
01622 z = l + len;
01623 mu = z + len;
01624 alpha = mu + len;
01625
01626 for (i = 0; i < len; i++) {
01627 a[i] = (double)y[i];
01628 }
01629
01630
01631 for (i = 0; i < end; i++) {
01632 h[i] = (double)x[i + 1] - (double)x[i];
01633 if (h[i] < 0.0) {
01634 cpl_free(a) ;
01635 return -1;
01636 }
01637 }
01638
01639
01640 for (n = 0, i = 1; i < end; i++, n++) {
01641
01642 alpha[i] = 3.0 * ((a[i+1] / h[i]) - (a[i] / h[n]) - (a[i] / h[i]) +
01643 (a[n] / h[n]));
01644 }
01645
01646
01647 l[0] = l[end] = 1.0;
01648 mu[0] = mu[end] = 0.0;
01649 z[0] = z[end] = 0.0;
01650 c[0] = c[end] = 0.0;
01651
01652
01653 for (n = 0, i = 1; i < end; i++, n++) {
01654
01655 l[i] = 2 * (h[i] + h[n]) - h[n] * mu[n];
01656 mu[i] = h[i] / l[i];
01657 z[i] = (alpha[i] - h[n] * z[n]) / l[i];
01658 }
01659 for (n = end, j = end - 1; j >= 0; j--, n--) {
01660
01661 c[j] = z[j] - mu[j] * c[n];
01662 b[j] = (a[n] - a[j]) / h[j] - h[j] * (c[n] + 2.0 * c[j]) / 3.0;
01663 d[j] = (c[n] - c[j]) / (3.0 * h[j]);
01664 }
01665
01666
01667 for (j = 0; j < splLen; j++) {
01668 v = (double)splX[j];
01669 splY[j] = (float)0;
01670
01671
01672 if ((v < (double)x[0]) || (v > (double)x[end])) {
01673 continue;
01674 }
01675
01676 loc = sinfo_function1d_search_value(x, len, (double)v, &found);
01677 if (found) {
01678 splY[j] = y[loc];
01679 } else {
01680 loc--;
01681 v -= (double)x[loc];
01682 splY[j] = (float)( a[loc] + v * (b[loc] + v * (c[loc] + v * d[loc])));
01683 }
01684 }
01685 cpl_free(a) ;
01686 return 0;
01687 }
01688
01689
01705
01706
01707 static int
01708 sinfo_function1d_search_value(
01709 double * x,
01710 int len,
01711 double key,
01712 int * foundPtr
01713 )
01714 {
01715 int high,
01716 low,
01717 middle;
01718
01719 low = 0;
01720 high = len - 1;
01721
01722 while (high >= low) {
01723 middle = (high + low) / 2;
01724 if (key > x[middle]) {
01725 low = middle + 1;
01726 } else if (key < x[middle]) {
01727 high = middle - 1;
01728 } else {
01729 *foundPtr = 1;
01730 return (middle);
01731 }
01732 }
01733 *foundPtr = 0;
01734 return (low);
01735 }
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772 static cpl_vector *
01773 sinfo_juha_vector_filter_median_create(
01774 const cpl_vector * v,
01775 int hw)
01776 {
01777 cpl_vector * filtered=NULL;
01778 double * row=NULL;
01779 int i, j, k, size;
01780 double temp;
01781
01782 size = cpl_vector_get_size(v);
01783 filtered = cpl_vector_new(size);
01784
01785 row = cpl_malloc((2*hw+1) * sizeof(double));
01786 for (i=0; i<size; i++) {
01787 k = 0;
01788 for (j=-hw; j<=hw; j++)
01789 if ( (i+j) >= 0 && (i+j) < size) {
01790 temp = cpl_vector_get (v, i+j);
01791 row[k] = temp;
01792 k++;
01793 }
01794 sinfo_tools_sort_double (row, k);
01795
01796 if (k%2 == 1)
01797 temp = row[k/2];
01798 else
01799 temp = row[k/2-1];
01800 cpl_vector_set (filtered, i, temp);
01801 }
01802 cpl_free(row);
01803 return filtered;
01804 }
01805
01806 #define CPL_PIX_STACK_SIZE 50
01807
01818
01819 static cpl_error_code sinfo_tools_sort_double(
01820 double * pix_arr,
01821 int n)
01822 {
01823 int i, ir, j, k, l;
01824 int * i_stack ;
01825 int j_stack ;
01826 double a ;
01827
01828
01829 cpl_ensure(pix_arr, CPL_ERROR_NULL_INPUT, CPL_ERROR_NULL_INPUT) ;
01830
01831 ir = n ;
01832 l = 1 ;
01833 j_stack = 0 ;
01834 i_stack = malloc(CPL_PIX_STACK_SIZE * sizeof(double)) ;
01835 for (;;) {
01836 if (ir-l < 7) {
01837 for (j=l+1 ; j<=ir ; j++) {
01838 a = pix_arr[j-1];
01839 for (i=j-1 ; i>=1 ; i--) {
01840 if (pix_arr[i-1] <= a) break;
01841 pix_arr[i] = pix_arr[i-1];
01842 }
01843 pix_arr[i] = a;
01844 }
01845 if (j_stack == 0) break;
01846 ir = i_stack[j_stack-- -1];
01847 l = i_stack[j_stack-- -1];
01848 } else {
01849 k = (l+ir) >> 1;
01850 SINFO_DOUBLE_SWAP(pix_arr[k-1], pix_arr[l])
01851 if (pix_arr[l] > pix_arr[ir-1]) {
01852 SINFO_DOUBLE_SWAP(pix_arr[l], pix_arr[ir-1])
01853 }
01854 if (pix_arr[l-1] > pix_arr[ir-1]) {
01855 SINFO_DOUBLE_SWAP(pix_arr[l-1], pix_arr[ir-1])
01856 }
01857 if (pix_arr[l] > pix_arr[l-1]) {
01858 SINFO_DOUBLE_SWAP(pix_arr[l], pix_arr[l-1])
01859 }
01860 i = l+1;
01861 j = ir;
01862 a = pix_arr[l-1];
01863 for (;;) {
01864 do i++; while (pix_arr[i-1] < a);
01865 do j--; while (pix_arr[j-1] > a);
01866 if (j < i) break;
01867 SINFO_DOUBLE_SWAP(pix_arr[i-1], pix_arr[j-1]);
01868 }
01869 pix_arr[l-1] = pix_arr[j-1];
01870 pix_arr[j-1] = a;
01871 j_stack += 2;
01872 if (j_stack > CPL_PIX_STACK_SIZE) {
01873
01874 free(i_stack);
01875 return CPL_ERROR_ILLEGAL_INPUT ;
01876 }
01877 if (ir-i+1 >= j-l) {
01878 i_stack[j_stack-1] = ir;
01879 i_stack[j_stack-2] = i;
01880 ir = j-1;
01881 } else {
01882 i_stack[j_stack-1] = j-1;
01883 i_stack[j_stack-2] = l;
01884 l = i;
01885 }
01886 }
01887 }
01888 free(i_stack) ;
01889 return CPL_ERROR_NONE ;
01890 }
01891
01892 static cpl_vector *
01893 sinfo_vector_filter_median_create(
01894 const cpl_vector * v,
01895 int hw)
01896 {
01897 cpl_vector * filtered;
01898 cpl_vector * row;
01899 int i, j, k, size;
01900 double temp;
01901
01902
01903 size = cpl_vector_get_size(v);
01904 filtered = cpl_vector_new(size);
01905
01906
01907 row = cpl_vector_new((2*hw+1));
01908 for (i=0; i<size; i++) {
01909 k = 0;
01910 for (j=-hw; j<=hw; j++)
01911 if ( (i+j) >= 0 && (i+j) < size) {
01912 temp = cpl_vector_get (v, i+j);
01913 cpl_vector_set(row,k,temp);
01914 k++;
01915 }
01916
01917
01918 cpl_vector_sort(row, +1);
01919 if (k%2 == 1) {
01920 temp = cpl_vector_get(row,k/2);
01921 }
01922 else {
01923 temp = cpl_vector_get(row,k/2-1);
01924 }
01925
01926 sinfo_msg("value = %g ", temp);
01927 cpl_vector_set (filtered, i, temp);
01928 }
01929 cpl_vector_delete(row);
01930 return filtered;
01931 }
01932
01933
01934
01935
01936 static double
01937 sinfo_image_get_median_window (const cpl_image *image,
01938 int llx, int lly, int urx, int ury)
01939 {
01940 cpl_image *window;
01941 float *data;
01942 double *array, median;
01943 int size, i,j;
01944
01945 window = cpl_image_extract (image, llx, lly, urx, ury);
01946 size = (urx-llx+1)*(ury-lly+1);
01947 data = cpl_image_get_data_float(window);
01948
01949 array = (double*)cpl_calloc ( size, sizeof(double));
01950 j = 0;
01951 for (i=0; i<size; i++)
01952 if (!isnan(data[i]))
01953 array[j++] = data[i];
01954
01955 if (j>0)
01956 sinfo_tools_sort_double (array, j);
01957
01958 if (j == 0 || 2*j<size)
01959 median = 0./0.;
01960 else if (j%2 == 1)
01961 median = array[j/2];
01962 else
01963 median = array[j/2-1];
01964
01965 cpl_image_delete (window);
01966 cpl_free (array);
01967
01968 return (median);
01969 }