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