00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "lamp_flats.h"
00023 #include "flat_ini_by_cpl.h"
00024 #include "sinfoni_pro_save.h"
00025
00026
00027
00028
00029
00030 int qc_get_cnt(cpl_frameset* on_set, cpl_frameset* of_set, flat_config* cfg);
00031 int lamp_flats_det_ncounts(cpl_frameset* raw, flat_config* config);
00032
00033 static struct {
00034
00035 double avg_on;
00036 double std_on;
00037 double avg_of;
00038 double std_of;
00039 double avg_di;
00040 double std_di;
00041 double nsat_on;
00042 double noise_on;
00043 double noise_of;
00044 double flux_on;
00045 double nsat;
00046
00047 } qc_lampflat;
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073 int lamp_flats (const char* plugin_id, cpl_parameterlist* config, cpl_frameset* sof)
00074 {
00075 const char* _id = "lamp_flats";
00076 flat_config * cfg =NULL;
00077 OneImage ** list_object ;
00078 OneImage ** list_dither_object ;
00079
00080 OneImage ** list_sky ;
00081 OneImage ** list_dither_sky;
00082
00083
00084
00085
00086
00087 OneImage ** im ;
00088 OneImage * norm_dith =NULL;
00089
00090
00091 OneImage * im_obj =NULL;
00092 OneImage * int_im =NULL;
00093 OneImage * int_im_dith =NULL;
00094
00095 OneImage * im_sky =NULL;
00096 OneImage * im_dither =NULL;
00097 OneImage * im_obj_sub =NULL;
00098 OneImage * im_dither_sub =NULL;
00099 OneImage * im_dither_sky =NULL;
00100 OneImage ** imMed ;
00101
00102
00103 OneImage * colImage =NULL;
00104 OneImage * mask_im =NULL;
00105 OneImage * norm =NULL;
00106
00107 OneImage * compImage =NULL;
00108 OneImage * maskImage =NULL;
00109 OneImage * threshIm =NULL;
00110 OneCube * cube_object =NULL;
00111 OneCube * cube_dither_object =NULL;
00112
00113 OneCube * cube_sky =NULL;
00114 OneCube * cube_dither_sky =NULL;
00115 char name[FILE_NAME_SZ];
00116
00117 int typ;
00118 Stats * stats =NULL;
00119 int i = 0;
00120 int* n=NULL;
00121 int nob =0;
00122 int nsky = 0;
00123 int nobjdith = 0;
00124 int nskydith = 0;
00125 int n_badpixels =0;
00126
00127 int pos =0;
00128
00129 float** slit_edges;
00130 float clean_mean =0.;
00131 float clean_stdev =0.;
00132 float mean_factor =0.;
00133 float val_x=0;
00134 float val_y=0;
00135
00136 char* outNameDither=NULL;
00137 char* name_list=NULL;
00138
00139 char* tbl_slitpos_name=NULL;
00140 cpl_table* tbl_slitpos=NULL;
00141 int* status=NULL;
00142
00143 int n_im_med=0;
00144 cpl_image* bp_img=NULL;
00145 cpl_image* obj_img=NULL;
00146 cpl_image* obj_imgw=NULL;
00147 cpl_image* dit_img=NULL;
00148
00149 cpl_frameset* raw=NULL;
00150
00151 cpl_table* qclog_tbl=NULL;
00152 char* key_value=NULL;
00153 int naxis1=0;
00154 int naxis2=0;
00155 double fpn_stdev1=0;
00156 double fpn_stdev2=0;
00157
00158
00159
00160
00161
00162
00163
00164
00165 raw=cpl_frameset_new();
00166 cfg = parse_cpl_input_flat(config,sof,&raw);
00167
00168 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00169 cpl_msg_error(_id,"error parsing cpl input");
00170 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00171 flat_free(cfg);
00172 cpl_frameset_delete(raw);
00173 return -1;
00174 }
00175
00176
00177
00178
00179 if (cfg == NULL)
00180 {
00181 cpl_msg_error (_id,"could not parse cpl input!\n") ;
00182 cpl_frameset_delete(raw);
00183 return -1 ;
00184 }
00185
00186 if (cfg->interpolInd == 1) {
00187 if(is_fits_file(cfg->mask) != 1) {
00188 cpl_msg_error(_id,"Input file %s is not FITS",cfg->mask);
00189 cpl_frameset_delete(raw);
00190 flat_free(cfg);
00191 return -1;
00192 }
00193
00194
00195 if (is_fits_file(cfg->slitposList) != 1) {
00196 cpl_msg_error(_id,"Input file %s is not FITS",cfg->slitposList);
00197 cpl_frameset_delete(raw);
00198 flat_free(cfg);
00199 return -1;
00200 }
00201 }
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211 cpl_msg_info(_id,"Takes clean mean of several images");
00212
00213 list_object = (OneImage**) cpl_calloc (cfg->nobj, sizeof(OneImage*)) ;
00214 if (list_object == NULL) {
00215 cpl_msg_error(_id,"could not allocate memory\n");
00216 cpl_frameset_delete(raw);
00217 flat_free(cfg);
00218 return -1;
00219 }
00220
00221
00222 if (cfg->contains_dither == 1) {
00223 list_dither_object=(OneImage**) cpl_calloc(cfg->nditherobj, sizeof(OneImage*));
00224 if (list_dither_object == NULL) {
00225 cpl_msg_error(_id,"could not allocate memory");
00226 destroy_list(list_object);
00227 cpl_frameset_delete(raw);
00228 flat_free(cfg);
00229 return -1;
00230 }
00231 }
00232
00233 if (cfg->contains_sky == 1) {
00234 list_sky = (OneImage**) cpl_calloc (cfg->noff, sizeof(OneImage*)) ;
00235 if (list_sky == NULL) {
00236 cpl_msg_error(_id,"could not allocate memory");
00237 destroy_list(list_object);
00238 cpl_frameset_delete(raw);
00239 flat_free(cfg);
00240 return -1;
00241 }
00242 }
00243
00244
00245
00246
00247 if (cfg->contains_dither == 1 && cfg->nditheroff > 0) {
00248 list_dither_sky=(OneImage**) cpl_calloc (cfg->nditheroff, sizeof(OneImage*));
00249 if (list_dither_sky == NULL) {
00250 cpl_msg_error(_id," could not allocate memory");
00251 destroy_list(list_sky);
00252 destroy_list(list_object);
00253 cpl_frameset_delete(raw);
00254 flat_free(cfg);
00255 return -1;
00256 }
00257 }
00258
00259
00260
00261
00262
00263 if (cfg->contains_dither == 0 && cfg->nditheroff > 0){
00264 cpl_msg_error(_id,"please use non-dithered off-frames, remove the 2!");
00265 destroy_list(list_sky);
00266 destroy_list(list_object);
00267 cpl_frameset_delete(raw);
00268 flat_free(cfg);
00269 return -1;
00270 }
00271
00272
00273
00274
00275
00276 im = (OneImage**) cpl_calloc (cfg -> nframes, sizeof(OneImage*)) ;
00277
00278
00279
00280
00281
00282
00283 for (i=0; i< cfg->nframes; i++) {
00284 strcpy(name,cfg->framelist[i]);
00285 if(is_fits_file(name) != 1) {
00286 cpl_msg_error(_id,"PP Input file %s %d is not FITS",name,i);
00287 if(im) cpl_free(im);
00288 destroy_list(list_sky);
00289 destroy_list(list_object);
00290 cpl_frameset_delete(raw);
00291 flat_free(cfg);
00292 return -1;
00293 }
00294 im[i] = load_image( name );
00295
00296
00297 }
00298
00299
00300 for (i=0; i< cfg->nframes; i++) {
00301 typ = intarray_get_value( cfg->frametype, i );
00302 pos = intarray_get_value( cfg->frameposition, i );
00303 if (im[i] == NULL) {
00304 cpl_msg_error(_id,"could not load image\n");
00305 destroy_list(list_sky);
00306 destroy_list(list_object);
00307 cpl_frameset_delete(raw);
00308 flat_free(cfg);
00309 return -1;
00310 }
00311 if (pos == 2) {
00312 if (typ == 1) {
00313 insert_image( list_object, im[i], nob );
00314
00315 nob = nob + 1;
00316 } else {
00317 insert_image( list_sky, im[i], nsky );
00318 nsky = nsky + 1 ;
00319 }
00320 } else {
00321 if (typ == 1) {
00322 insert_image( list_dither_object, im[i], nobjdith );
00323
00324 nobjdith = nobjdith + 1;
00325 } else {
00326 insert_image( list_dither_sky, im[i], nskydith );
00327 nskydith = nskydith + 1 ;
00328 }
00329 }
00330
00331 }
00332
00333
00334
00335 if (nob != cfg->nobj || cfg->noff != nsky ||
00336 nobjdith != cfg->nditherobj || nskydith != cfg->nditheroff) {
00337 cpl_msg_error(_id,"something is wrong with the number of the different types of frames");
00338
00339 for (i=0; i< cfg->nframes; i++) {
00340 destroy_image(im[i]);
00341 }
00342 cpl_free(im);
00343 destroy_list(list_sky);
00344 destroy_list(list_object);
00345 cpl_frameset_delete(raw);
00346 flat_free(cfg);
00347 return -1;
00348 }
00349
00350
00351
00352
00353 cpl_msg_info(_id,"Creates and fills cubes with the different image lists");
00354 cube_object = list_make_cube(list_object, nob);
00355 if (cube_object == NULL) {
00356 cpl_msg_error(_id,"could not create data cube!");
00357 cpl_free(im);
00358 destroy_list(list_sky);
00359 destroy_list(list_object);
00360 cpl_frameset_delete(raw);
00361 flat_free(cfg);
00362 return -1;
00363 }
00364 cpl_free(list_object);
00365 if (cfg->contains_dither == 1) {
00366 cube_dither_object = list_make_cube(list_dither_object, nobjdith);
00367 if (cube_dither_object == NULL) {
00368 cpl_msg_error(_id,"could not create data cube!\n");
00369 cpl_free(im);
00370 destroy_list(list_sky);
00371 destroy_list(list_object);
00372 cpl_frameset_delete(raw);
00373 flat_free(cfg);
00374 return -1;
00375 }
00376 cpl_free(list_dither_object);
00377 }
00378 if (cfg->contains_sky == 1 && nsky > 0) {
00379 cube_sky = list_make_cube(list_sky, nsky);
00380 if (cube_sky == NULL) {
00381 cpl_msg_error(_id,"could not create data cube!\n");
00382 cpl_free(im);
00383 destroy_list(list_sky);
00384 destroy_list(list_object);
00385 cpl_frameset_delete(raw);
00386 flat_free(cfg);
00387 return -1;
00388 }
00389 cpl_free(list_sky);
00390 }
00391 if (cfg->contains_dither == 1 && nskydith > 0) {
00392 cube_dither_sky = list_make_cube(list_dither_sky, nskydith);
00393 if (cube_dither_sky == NULL) {
00394 cpl_msg_error(_id,"could not create data cube!\n");
00395 cpl_free(im);
00396 destroy_cube(cube_sky);
00397 destroy_cube(cube_object);
00398 cpl_frameset_delete(raw);
00399 flat_free(cfg);
00400 return -1;
00401 }
00402 cpl_free(list_dither_sky);
00403 }
00404
00405
00406
00407 cpl_msg_info(_id,"Takes the average of the different cubes");
00408 if (cfg->loReject*cfg->nobj < 1. && cfg->hiReject *cfg->nobj < 1.) {
00409 im_obj = averageCubeToImage( cube_object );
00410 if (im_obj == NULL) {
00411 cpl_msg_error(_id,"averageCubeToImage failed\n" );
00412 return -1;
00413 }
00414 } else {
00415 im_obj = average_with_rejection( cube_object,
00416 cfg->loReject, cfg->hiReject );
00417 if (im_obj == NULL) {
00418 cpl_msg_error(_id,"average_with_rejection failed\n" );
00419 return -1;
00420 }
00421 }
00422 destroy_cube(cube_object);
00423
00424 if (cfg->contains_sky == 1) {
00425 if (cfg->loReject * nsky < 1. && cfg->hiReject * nsky < 1.) {
00426 im_sky = averageCubeToImage( cube_sky );
00427 if (im_sky == NULL) {
00428 cpl_msg_error(_id, "averageCubeToImage failed\n" );
00429 cpl_free(im);
00430 destroy_list(list_sky);
00431 destroy_list(list_object);
00432 cpl_frameset_delete(raw);
00433 flat_free(cfg);
00434 return -1;
00435 }
00436 } else {
00437 im_sky = average_with_rejection( cube_sky,
00438 cfg->loReject, cfg->hiReject );
00439
00440 if (im_sky == NULL) {
00441 cpl_msg_error(_id, "average_with_rejection failed\n" );
00442 cpl_free(im);
00443 destroy_list(list_sky);
00444 destroy_list(list_object);
00445 cpl_frameset_delete(raw);
00446 flat_free(cfg);
00447 return -1;
00448 }
00449 }
00450 destroy_cube(cube_sky);
00451 }
00452
00453
00454
00455
00456 if (cfg->contains_dither == 1) {
00457 if (cfg->loReject*nobjdith < 1. && cfg->hiReject * nobjdith < 1.) {
00458 im_dither = averageCubeToImage( cube_dither_object );
00459 if (im_dither == NULL) {
00460 cpl_msg_error (_id, "averageCubeToImage failed\n" );
00461 cpl_free(im);
00462 destroy_list(list_sky);
00463 destroy_list(list_object);
00464 cpl_frameset_delete(raw);
00465 flat_free(cfg);
00466 return -1;
00467 }
00468 } else {
00469 im_dither = average_with_rejection( cube_dither_object,
00470 cfg->loReject, cfg->hiReject );
00471
00472 if (im_dither == NULL) {
00473 cpl_msg_error(_id, "average_with_rejection failed\n" );
00474 cpl_free(im);
00475 destroy_list(list_sky);
00476 destroy_list(list_object);
00477 cpl_frameset_delete(raw);
00478 flat_free(cfg);
00479 return -1;
00480 }
00481 }
00482 destroy_cube(cube_dither_object);
00483 }
00484
00485
00486 if (cfg->contains_dither == 1 && nskydith > 0 ) {
00487 if (cfg->loReject*nskydith < 1. && cfg->hiReject*nskydith < 1.) {
00488 im_dither_sky = averageCubeToImage( cube_dither_sky );
00489 if (im_dither_sky == NULL) {
00490 cpl_msg_error(_id, "averageCubeToImage failed\n" );
00491 cpl_free(im);
00492 destroy_list(list_sky);
00493 destroy_list(list_object);
00494 cpl_frameset_delete(raw);
00495 flat_free(cfg);
00496 return -1;
00497 }
00498 } else {
00499 im_dither_sky = average_with_rejection( cube_dither_sky,
00500 cfg->loReject, cfg->hiReject );
00501
00502 if (im_dither_sky == NULL) {
00503 cpl_msg_error(_id, "average_with_rejection failed\n" );
00504 cpl_free(im);
00505 destroy_list(list_sky);
00506 destroy_list(list_object);
00507 cpl_frameset_delete(raw);
00508 flat_free(cfg);
00509 return -1;
00510 }
00511 }
00512 destroy_cube(cube_dither_sky);
00513
00514 }
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527 cpl_msg_info(_id,"Subtracts the resulting off-frame (sky) from the on-frame");
00528 if (cfg->contains_sky == 1) {
00529 im_obj_sub = sub_image(im_obj, im_sky);
00530 if (im_obj_sub == NULL) {
00531 cpl_msg_error(_id,"could not sub_image\n");
00532 cpl_free(im);
00533 destroy_list(list_sky);
00534 destroy_list(list_object);
00535 cpl_frameset_delete(raw);
00536 flat_free(cfg);
00537 return -1;
00538 }
00539 destroy_image(im_obj);
00540 if (((cfg->contains_dither == 1) && (nskydith > 0)) ||
00541 (cfg->contains_dither == 0)) {
00542 destroy_image(im_sky);
00543 }
00544 im_obj = im_obj_sub;
00545 }
00546 if (cfg->contains_dither == 1 && nskydith > 0) {
00547 im_dither_sub = sub_image(im_dither, im_dither_sky);
00548 if (im_dither_sub == NULL) {
00549 cpl_msg_error(_id,"could not sub_image\n");
00550 cpl_free(im);
00551 destroy_list(list_sky);
00552 destroy_list(list_object);
00553 cpl_frameset_delete(raw);
00554 flat_free(cfg);
00555 return -1;
00556 }
00557 destroy_image(im_dither);
00558 destroy_image(im_dither_sky);
00559 im_dither = im_dither_sub;
00560 } else if (cfg->contains_dither == 1 &&
00561 nskydith == 0 &&
00562 cfg->contains_sky == 1) {
00563 im_dither_sub = sub_image(im_dither, im_sky);
00564 if (im_dither_sub == NULL) {
00565 cpl_msg_error(_id,"could not sub_image\n");
00566 cpl_free(im);
00567 destroy_list(list_sky);
00568 destroy_list(list_object);
00569 cpl_frameset_delete(raw);
00570 flat_free(cfg);
00571 return -1;
00572 }
00573 destroy_image(im_dither);
00574 destroy_image(im_sky);
00575 im_dither = im_dither_sub;
00576 }
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586 cpl_msg_info(_id,"Generating a static bad pixel mask");
00587 n_im_med = cfg->iterations+1;
00588
00589 imMed=(OneImage**) cpl_calloc(n_im_med, sizeof(OneImage*));
00590
00591 if (cfg->badInd == 1) {
00592 cpl_msg_info(_id,"removes the intensity tilt from every column and");
00593 cpl_msg_info(_id,"computes the standard deviation on a rectangular zone");
00594
00595
00596 colImage = colTilt( im_obj, cfg->sigmaFactor );
00597
00598 if (colImage == NULL) {
00599 cpl_msg_error(_id, "colTilt failed\n" );
00600 return -1;
00601 }
00602
00603 stats = imageStatsOnRectangle(colImage,
00604 cfg->badLoReject,
00605 cfg->badHiReject,
00606 cfg->llx,
00607 cfg->lly,
00608 cfg->urx,
00609 cfg->ury);
00610
00611
00612 if (stats == NULL) {
00613 cpl_msg_error (_id,"get_image_stats_on_vig failed\n");
00614 return -1;
00615 }
00616 clean_mean = Stats_get_cleanmean(stats);
00617 clean_stdev = Stats_get_cleanstdev(stats);
00618
00619
00620 if (cfg->threshInd == 1) {
00621 threshIm = threshImage(colImage,
00622 clean_mean-mean_factor*clean_stdev,
00623 clean_mean+mean_factor*clean_stdev);
00624 if (threshIm == NULL) {
00625 cpl_msg_error (_id, " threshImage failed\n" );
00626 return -1;
00627 }
00628 }
00629
00630
00631 if (cfg->threshInd == 0) {
00632 threshIm = colImage;
00633 }
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643 if ((imMed[0]= medianImage(threshIm, -cfg->factor*clean_stdev)) == NULL) {
00644 cpl_msg_error(_id, " medianImage failed\n" );
00645 return -1;
00646 }
00647
00648
00649
00650 for (i=1; i< cfg->iterations+1; i++) {
00651 if ((imMed[i]=medianImage(imMed[i-1],
00652 -cfg->factor*clean_stdev)) == NULL) {
00653 cpl_msg_error(_id, "medianImage failed\n" );
00654 return -1;
00655 }
00656 }
00657
00658
00659
00660
00661
00662
00663 compImage = compareImages(threshIm, imMed[cfg->iterations], im_obj);
00664 if (compImage == NULL) {
00665 cpl_msg_error(_id, "compareImages failed\n" );
00666
00667 for (i=0; i< cfg->iterations+1; i++) {
00668 destroy_image(imMed[i]);
00669 }
00670 delStats(stats);
00671 destroy_image(im_obj);
00672 cpl_free(imMed);
00673 destroy_image(colImage);
00674 cpl_free(im);
00675 cpl_frameset_delete(raw);
00676 flat_free(cfg);
00677 return -1;
00678 }
00679
00680
00681
00682
00683
00684 n = new_int_array(1);
00685 maskImage = promoteImageToMask( compImage, n );
00686 if (maskImage == NULL) {
00687 cpl_msg_error( _id,"error in promoteImageToMask" );
00688 destroy_image(compImage);
00689 for (i=0; i< cfg->iterations+1; i++) {
00690 destroy_image(imMed[i]);
00691 }
00692 delStats(stats);
00693 destroy_image(im_obj);
00694 cpl_free(imMed);
00695 destroy_image(colImage);
00696 cpl_free(im);
00697 cpl_frameset_delete(raw);
00698 flat_free(cfg);
00699 return -1;
00700 }
00701
00702
00703 n_badpixels = intarray_get_value(n, 0);
00704 cpl_msg_info(_id,"No of bad pixels: %d", n_badpixels);
00705 bp_img=cpl_image_wrap_float(maskImage->lx,maskImage->ly,maskImage->data);
00706
00707
00708
00709 qclog_tbl = sinfoni_qclog_init(1);
00710 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00711 sprintf(key_value,"%d",n_badpixels);
00712 sinfoni_qclog_add(qclog_tbl,0,"QC BP_MAP NBADPIX","CPL_TYPE_INT",
00713 key_value,"No of bad pixels");
00714
00715 if(-1 == sinfoni_pro_save_ima(bp_img,raw,sof,cfg->maskname,
00716 PRO_BP_MAP,qclog_tbl,plugin_id,config)) {
00717 cpl_msg_error(_id,"cannot save ima %s", cfg->maskname);
00718
00719 cpl_table_delete(qclog_tbl);
00720
00721 for (i=0; i< n_im_med; i++) {
00722 destroy_image(imMed[i]);
00723 }
00724 cpl_free(imMed);
00725 delStats(stats);
00726 destroy_intarray(n);
00727 destroy_image(threshIm);
00728 if (cfg->threshInd == 1) {
00729 destroy_image(colImage);
00730 }
00731 destroy_image(compImage);
00732 destroy_image(maskImage);
00733
00734 if (im) {
00735 cpl_free(im);
00736 }
00737
00738 cpl_free(name_list);
00739 cpl_free(outNameDither);
00740 flat_free(cfg);
00741
00742 return -1;
00743 }
00744
00745 cpl_table_delete(qclog_tbl);
00746 cpl_free(key_value);
00747
00748
00749
00750
00751 delStats(stats);
00752 destroy_intarray(n);
00753 destroy_image(threshIm);
00754 if (cfg->threshInd == 1) {
00755 destroy_image(colImage);
00756 }
00757 destroy_image(compImage);
00758 destroy_image(maskImage);
00759
00760
00761 for (i=0; i< cfg->iterations+1; i++) {
00762 destroy_image(imMed[i]);
00763 }
00764
00765 cpl_image_unwrap(bp_img);
00766
00767
00768 }
00769 cpl_free(imMed);
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781 cpl_msg_info(_id,"Creates a Master flat field");
00782 if (cfg->interpolInd == 1) {
00783 mask_im = load_image(cfg->mask);
00784 if (mask_im == NULL) {
00785 cpl_msg_error(_id, "could not load static bad pixel mask\n" );
00786 return -1;
00787 }
00788
00789
00790 slit_edges = new_2Dfloat_array(32, 2) ;
00791
00792
00793 if(is_fits_file(cfg->slitposList) !=1 ) {
00794 cpl_msg_error(_id,"Input file %s is not FITS", cfg->slitposList);
00795 return -1;
00796 }
00797 tbl_slitpos_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00798 strcpy(tbl_slitpos_name,cfg->slitposList);
00799 tbl_slitpos = cpl_table_load(tbl_slitpos_name,1,0);
00800 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00801 cpl_msg_error(_id,"error loading tbl %s",tbl_slitpos_name);
00802 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00803 return -1;
00804 }
00805
00806 for (i =0 ; i< 32; i++){
00807
00808 val_x=cpl_table_get_double(tbl_slitpos,"pos1",i,status);
00809 val_y=cpl_table_get_double(tbl_slitpos,"pos2",i,status);
00810 array2D_set_value(slit_edges,val_x,i,0);
00811 array2D_set_value(slit_edges,val_y,i,1);
00812 }
00813 cpl_table_delete(tbl_slitpos);
00814 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00815 cpl_msg_error(_id,"reading tbl %s",tbl_slitpos_name);
00816 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00817 return -1;
00818 }
00819 cpl_free(tbl_slitpos_name);
00820
00821 int_im = interpolSourceImage (im_obj, mask_im, cfg->maxRad, slit_edges);
00822 if (int_im == NULL) {
00823 cpl_msg_error(_id, "could not carry out interpolSourceImage\n" );
00824 return -1;
00825 }
00826 destroy_image(im_obj);
00827 norm = normalize_to_central_pixel(int_im);
00828 if (norm == NULL) {
00829 cpl_msg_error(_id, "could not normalize flatfield\n" );
00830 return -1;
00831 }
00832 destroy_image(int_im);
00833 im_obj = norm;
00834 if (cfg->contains_dither == 1) {
00835 int_im_dith = interpolSourceImage (im_dither, mask_im, cfg->maxRad,
00836 slit_edges);
00837 if (int_im_dith == NULL) {
00838 cpl_msg_error(_id, "could not carry out interpolSourceImage\n" );
00839 return -1;
00840 }
00841 destroy_image(im_dither);
00842 norm_dith = normalize_to_central_pixel(int_im_dith);
00843 if (norm_dith == NULL) {
00844 cpl_msg_error(_id, "could not normalize flatfield\n" );
00845 return -1;
00846 }
00847 destroy_image(int_im_dith);
00848 im_dither = norm_dith;
00849 }
00850 destroy_2Darray(slit_edges, 32);
00851 destroy_image(mask_im);
00852
00853 }
00854
00855 if (cfg->interpolInd != 1) {
00856 norm = normalize_to_central_pixel(im_obj);
00857 if (norm == NULL) {
00858 cpl_msg_error(_id, "could not normalize flatfield\n" );
00859 return -1;
00860 }
00861 destroy_image(im_obj);
00862 im_obj = norm;
00863
00864 if (cfg->contains_dither == 1) {
00865 norm_dith = normalize_to_central_pixel(im_dither);
00866 if (norm_dith == NULL) {
00867 cpl_msg_error(_id, "could not normalize flatfield\n" );
00868 return -1;
00869 }
00870 destroy_image(im_dither);
00871 im_dither = norm_dith;
00872 }
00873 }
00874
00875 obj_imgw=cpl_image_wrap_float(im_obj->lx,im_obj->ly,im_obj->data);
00876 obj_img=cpl_image_duplicate(obj_imgw);
00877 cpl_image_unwrap(obj_imgw);
00878
00879
00880 naxis1=cpl_image_get_size_x(obj_img);
00881 naxis2=cpl_image_get_size_y(obj_img);
00882
00883
00884 if(cfg->qc_fpn_xmin1 < 1) {
00885 cpl_msg_error(_id,"qc_ron_xmin < 1");
00886 return -1;
00887 }
00888
00889 if(cfg->qc_fpn_xmax1 > naxis1) {
00890 cpl_msg_error(_id,"qc_ron_xmax < %d",naxis1);
00891 return -1;
00892 }
00893
00894 if(cfg->qc_fpn_ymin1 < 1) {
00895 cpl_msg_error(_id,"qc_ron_ymin < 1");
00896 return -1;
00897 }
00898
00899 if(cfg->qc_fpn_ymax1 > naxis2) {
00900 cpl_msg_error(_id,"qc_ron_ymax < %d",naxis2);
00901 return -1;
00902 }
00903
00904
00905
00906
00907 fpn_stdev1 = cpl_image_get_stdev_window(obj_img,
00908 cfg->qc_fpn_xmin1,
00909 cfg->qc_fpn_ymin1,
00910 cfg->qc_fpn_xmax1,
00911 cfg->qc_fpn_ymax1);
00912
00913
00914 if(cfg->qc_fpn_xmin2 < 1) {
00915 cpl_msg_error(_id,"qc_ron_xmin < %d",1);
00916 return -1;
00917 }
00918
00919
00920 if(cfg->qc_fpn_xmax2 > naxis1) {
00921 cpl_msg_error(_id,"qc_ron_xmax < %d",naxis1);
00922 return -1;
00923 }
00924
00925 if(cfg->qc_fpn_ymin2 < 1) {
00926 cpl_msg_error(_id,"qc_ron_ymin < 1");
00927 return -1;
00928 }
00929
00930 if(cfg->qc_fpn_ymax2 > naxis2) {
00931 cpl_msg_error(_id,"qc_ron_ymax < %d",naxis2);
00932 return -1;
00933 }
00934
00935
00936
00937 fpn_stdev2 = cpl_image_get_stdev_window(obj_img,
00938 cfg->qc_fpn_xmin2,
00939 cfg->qc_fpn_ymin2,
00940 cfg->qc_fpn_xmax2,
00941 cfg->qc_fpn_ymax2);
00942
00943
00944 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00945 cpl_msg_error(_id,"error with stdev_subw ");
00946 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00947 return -1;
00948 }
00949
00950
00951
00952
00953
00954 lamp_flats_det_ncounts(raw,cfg);
00955
00956
00957
00958
00959
00960 qclog_tbl = sinfoni_qclog_init(5);
00961 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00962 sprintf(key_value,"%g",qc_lampflat.avg_di);
00963 sinfoni_qclog_add(qclog_tbl,0,"QC SPECFLAT NCNTSAVG","CPL_TYPE_DOUBLE",
00964 key_value,"Average counts");
00965
00966 sprintf(key_value,"%g",qc_lampflat.std_di);
00967 sinfoni_qclog_add(qclog_tbl,1,"QC SPECFLAT NCNTSSTD","CPL_TYPE_DOUBLE",
00968 key_value,"Stdev counts");
00969
00970 sprintf(key_value,"%g",qc_lampflat.avg_of);
00971 sinfoni_qclog_add(qclog_tbl,2,"QC SPECFLAT OFFFLUX","CPL_TYPE_DOUBLE",
00972 key_value,"Average flux off frames");
00973
00974 sprintf(key_value,"%f",fpn_stdev1);
00975 sinfoni_qclog_add(qclog_tbl,3,"QC LFLAT FPN1","CPL_TYPE_DOUBLE",
00976 key_value,"Fixed Patter Noise of combined frames");
00977
00978 sprintf(key_value,"%f",fpn_stdev2);
00979 sinfoni_qclog_add(qclog_tbl,4,"QC LFLAT FPN2","CPL_TYPE_DOUBLE",
00980 key_value,"Fixed Patter Noise of combined frames");
00981
00982
00983
00984
00985
00986
00987 if(-1 == sinfoni_pro_save_ima(obj_img,raw,sof,cfg->outName,
00988 PRO_MASTER_FLAT_LAMP,qclog_tbl,plugin_id,config)) {
00989 cpl_msg_error(_id,"cannot save ima %s", cfg->outName);
00990 cpl_table_delete(qclog_tbl);
00991 destroy_image(im_obj);
00992 cpl_image_delete(obj_img);
00993 cpl_free(key_value);
00994 if (im) {
00995 cpl_free(im);
00996 }
00997 flat_free(cfg);
00998 cpl_frameset_delete(raw);
00999 return -1;
01000 }
01001 cpl_table_delete(qclog_tbl);
01002 destroy_image(im_obj);
01003 cpl_free(key_value);
01004
01005
01006
01007
01008 name_list = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char));
01009 outNameDither = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char));
01010
01011 if (cfg->contains_dither == 1) {
01012
01013 if (strstr(cfg->outName, ".fits" ) != NULL ) {
01014
01015 sprintf(name_list, "%s%s", get_rootname(cfg->outName), "_dither");
01016 strcpy(outNameDither,name_list);
01017 strcat(outNameDither,strstr(cfg->outName,".fits"));
01018
01019 } else {
01020 strcpy(outNameDither,cfg->outName);
01021 strcat(outNameDither,"_dither");
01022 }
01023
01024
01025
01026 dit_img=cpl_image_wrap_float(im_dither->lx,
01027 im_dither->ly,
01028 im_dither->data);
01029
01030
01031 naxis1=cpl_image_get_size_x(dit_img);
01032 naxis2=cpl_image_get_size_y(dit_img);
01033
01034
01035 if(cfg->qc_fpn_xmin1 < 1) {
01036 cpl_msg_error(_id,"qc_ron_xmin1 < 1");
01037 return -1;
01038 }
01039
01040 if(cfg->qc_fpn_xmax1 > naxis1) {
01041 cpl_msg_error(_id,"qc_ron_xmax1 < %d",naxis1);
01042 return -1;
01043 }
01044
01045 if(cfg->qc_fpn_ymin1 < 1) {
01046 cpl_msg_error(_id,"qc_ron_ymin1 < 1");
01047 return -1;
01048 }
01049
01050 if(cfg->qc_fpn_ymax1 > naxis2) {
01051 cpl_msg_error(_id,"qc_ron_ymax1 < %d",naxis2);
01052 return -1;
01053 }
01054
01055
01056 fpn_stdev1 = cpl_image_get_stdev_window(dit_img,
01057 cfg->qc_fpn_xmin1,
01058 cfg->qc_fpn_ymin1,
01059 cfg->qc_fpn_xmax1,
01060 cfg->qc_fpn_ymax1);
01061
01062
01063
01064
01065
01066
01067
01068 if(cfg->qc_fpn_xmin2 < 1) {
01069 cpl_msg_error(_id,"qc_ron_xmin2 < 1");
01070 return -1;
01071 }
01072
01073 if(cfg->qc_fpn_xmax2 > naxis1) {
01074 cpl_msg_error(_id,"qc_ron_xmax2 < %d",naxis1);
01075 return -1;
01076 }
01077
01078 if(cfg->qc_fpn_ymin2 < 1) {
01079 cpl_msg_error(_id,"qc_ron_ymin2 < 1");
01080 return -1;
01081 }
01082
01083 if(cfg->qc_fpn_ymax2 > naxis2) {
01084 cpl_msg_error(_id,"qc_ron_ymax2 < %d",naxis2);
01085 return -1;
01086 }
01087
01088
01089 fpn_stdev2 = cpl_image_get_stdev_window(dit_img,
01090 cfg->qc_fpn_xmin2,
01091 cfg->qc_fpn_ymin2,
01092 cfg->qc_fpn_xmax2,
01093 cfg->qc_fpn_ymax2);
01094
01095
01096 lamp_flats_det_ncounts(raw,cfg);
01097 qclog_tbl = sinfoni_qclog_init(5);
01098
01099 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
01100 sprintf(key_value,"%g",qc_lampflat.avg_di);
01101
01102 sinfoni_qclog_add(qclog_tbl,0,"QC SPECFLAT NCNTSAVG","CPL_TYPE_DOUBLE",
01103 key_value,"Average counts");
01104
01105 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
01106 sprintf(key_value,"%g",qc_lampflat.std_di);
01107 sinfoni_qclog_add(qclog_tbl,1,"QC SPECFLAT NCNTSSTD","CPL_TYPE_DOUBLE",
01108 key_value,"Stdev counts");
01109
01110 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
01111 sprintf(key_value,"%g",qc_lampflat.avg_of);
01112
01113 sinfoni_qclog_add(qclog_tbl,2,"QC SPECFLAT OFFFLUX","CPL_TYPE_DOUBLE",
01114 key_value,"Average flux off frames");
01115
01116 sprintf(key_value,"%f",fpn_stdev1);
01117 sinfoni_qclog_add(qclog_tbl,3,"QC LFLAT FPN1","CPL_TYPE_DOUBLE",
01118 key_value,"Fixed Patter Noise of combined frames");
01119
01120 sprintf(key_value,"%f",fpn_stdev2);
01121 sinfoni_qclog_add(qclog_tbl,4,"QC LFLAT FPN2","CPL_TYPE_DOUBLE",
01122 key_value,"Fixed Patter Noise of combined frames");
01123
01124 if(-1 == sinfoni_pro_save_ima(dit_img,raw,sof,outNameDither,
01125 PRO_MASTER_FLAT_LAMP,qclog_tbl,plugin_id,config)) {
01126 cpl_msg_error(_id,"cannot save ima %s", outNameDither);
01127
01128 cpl_free(imMed);
01129
01130
01131
01132
01133
01134
01135
01136
01137 if (im) {
01138 cpl_free(im);
01139 }
01140
01141 cpl_free(name_list);
01142 cpl_free(outNameDither);
01143 flat_free(cfg);
01144
01145 return -1;
01146 }
01147 cpl_table_delete(qclog_tbl);
01148 destroy_image(im_dither);
01149
01150 }
01151
01152
01153
01154
01155
01156
01157 if(name_list) cpl_free(name_list);
01158 if(outNameDither) cpl_free(outNameDither);
01159 cpl_image_delete(obj_img);
01160 cpl_free(im);
01161 cpl_frameset_delete(raw);
01162 flat_free(cfg);
01163 return 0;
01164
01165
01166
01167 }
01168
01169 int
01170 lamp_flats_det_ncounts(cpl_frameset* raw, flat_config* cfg)
01171 {
01172 int i=0;
01173 int j=0;
01174
01175 int nraw=0;
01176 int non=0;
01177 int noff=0;
01178
01179 double mjd_on=0;
01180 double mjd_of=0;
01181 double mjd_of_frm=0;
01182
01183 const char* _id = "lamp_flats_det_ncounts";
01184 char filename[FILE_NAME_SZ];
01185
01186
01187 cpl_frame* frm=NULL;
01188 cpl_frame* frm_dup=NULL;
01189 cpl_frame* on_frm=NULL;
01190 cpl_frame* of_frm=NULL;
01191 cpl_frame* tmp_of_frm=NULL;
01192
01193
01194 cpl_frameset* on_set=NULL;
01195 cpl_frameset* of_set=NULL;
01196 cpl_frameset* wrk_set=NULL;
01197
01198 on_set=cpl_frameset_new();
01199 of_set=cpl_frameset_new();
01200
01201 nraw = cpl_frameset_get_size(raw);
01202
01203 for (i=0; i< nraw; i++) {
01204 frm = cpl_frameset_get_frame(raw,i);
01205 frm_dup = cpl_frame_duplicate(frm);
01206 if(sinfoni_frame_is_on(frm) == 1) {
01207 cpl_frameset_insert(on_set,frm_dup);
01208 non++;
01209 } else {
01210 cpl_frameset_insert(of_set,frm_dup);
01211 noff++;
01212 }
01213 }
01214
01215
01216 if (non == noff) {
01217 qc_get_cnt(on_set,of_set,cfg);
01218
01219 } else if (non == 0) {
01220 cpl_msg_info(_id,"non == 0");
01221 cpl_msg_warning(_id,"QC SPECFLAT NCNTAVG/NCTNTSTD/OFFFLUX=0 ");
01222
01223 } else if ( noff == 0 ) {
01224 cpl_msg_info(_id,"noff == 0");
01225 cpl_msg_warning(_id,"QC SPECFLAT NCNTAVG/NCTNTSTD/OFFFLUX=0 ");
01226
01227 } else {
01228
01229 cpl_msg_warning(_id,"non != noff, => QC SPECFLAT NCNTAVG/NCTNTSTD/OFFFLUX=0 ");
01230
01231 for (i=0;i<non;i++) {
01232 wrk_set=cpl_frameset_new();
01233 on_frm=cpl_frameset_get_frame(on_set,i);
01234 mjd_on=sinfoni_get_mjd_obs(on_frm);
01235 of_frm=cpl_frameset_get_frame(of_set,0);
01236 mjd_of=sinfoni_get_mjd_obs(of_frm);
01237 strcpy(filename,cpl_frame_get_filename(of_frm));
01238 for (j=1;j<noff;j++) {
01239 tmp_of_frm = cpl_frameset_get_frame(of_set,j);
01240 mjd_of_frm = sinfoni_get_mjd_obs(tmp_of_frm);
01241
01242 if(1000.*(mjd_of_frm-mjd_on)*(mjd_of_frm-mjd_on) <
01243 1000.*(mjd_of- mjd_on)*(mjd_of- mjd_on) ) {
01244 mjd_of=mjd_of_frm;
01245 of_frm=cpl_frame_duplicate(tmp_of_frm);
01246 }
01247 }
01248 strcpy(filename,cpl_frame_get_filename(of_frm));
01249 frm_dup=cpl_frame_duplicate(of_frm);
01250 cpl_frameset_insert(wrk_set,frm_dup);
01251 strcpy(filename,cpl_frame_get_filename(of_frm));
01252 }
01253
01254 qc_get_cnt(on_set,wrk_set,cfg);
01255
01256 }
01257
01258 cpl_frameset_delete(wrk_set);
01259 cpl_frameset_delete(on_set);
01260 cpl_frameset_delete(of_set);
01261 return 0;
01262
01263
01264 }
01265
01266 int
01267 qc_get_cnt(cpl_frameset* on_set, cpl_frameset* of_set, flat_config* cfg)
01268 {
01269
01270 int i=0;
01271 int nsat=0;
01272 int non=0;
01273 int nof=0;
01274 int nfr=0;
01275
01276 char name[FILE_NAME_SZ];
01277 cpl_vector* vec_on=NULL;
01278 cpl_vector* vec_of=NULL;
01279 cpl_vector* vec_di=NULL;
01280 cpl_vector* vec_nsat=NULL;
01281 cpl_frame* on_frm=NULL;
01282 cpl_frame* of_frm=NULL;
01283
01284 cpl_image* dif_ima=NULL;
01285 cpl_image* on_ima=NULL;
01286 cpl_image* of_ima=NULL;
01287 cpl_image* tmp_ima=NULL;
01288
01289 double med=0;
01290 non = cpl_frameset_get_size(on_set);
01291 nof = cpl_frameset_get_size(of_set);
01292 nfr = (non <= nof) ? non : nof;
01293 vec_on = cpl_vector_new(nfr);
01294 vec_of = cpl_vector_new(nfr);
01295 vec_di = cpl_vector_new(nfr);
01296 vec_nsat = cpl_vector_new(nfr);
01297
01298
01299 for (i=0; i< nfr; i++) {
01300 on_frm = cpl_frameset_get_frame(on_set,i);
01301 strcpy(name,cpl_frame_get_filename(on_frm));
01302 on_ima = cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
01303 med= cpl_image_get_median(on_ima);
01304 cpl_vector_set(vec_on,i,med);
01305
01306 tmp_ima = cpl_image_duplicate(on_ima);
01307 cpl_image_threshold(tmp_ima,CPL_PIXEL_MINVAL,
01308 cfg->qc_thresh_max,0,1);
01309 nsat=cpl_image_get_flux(tmp_ima);
01310 cpl_vector_set(vec_nsat,i,nsat);
01311
01312
01313 of_frm = cpl_frameset_get_frame(of_set,i);
01314 strcpy(name,cpl_frame_get_filename(of_frm));
01315 of_ima = cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
01316 med= cpl_image_get_median(of_ima);
01317 cpl_vector_set(vec_of,i,med);
01318 dif_ima = cpl_image_subtract_create(on_ima,of_ima);
01319 med= cpl_image_get_median(dif_ima);
01320 cpl_vector_set(vec_di,i,med);
01321
01322 cpl_image_delete(on_ima);
01323 cpl_image_delete(of_ima);
01324 cpl_image_delete(dif_ima);
01325 cpl_image_delete(tmp_ima);
01326 }
01327 qc_lampflat.avg_on=cpl_vector_get_mean(vec_on);
01328 qc_lampflat.avg_of=cpl_vector_get_mean(vec_of);
01329 qc_lampflat.avg_di=cpl_vector_get_mean(vec_di);
01330 if(nfr > 1 ) {
01331 qc_lampflat.std_on=cpl_vector_get_stdev(vec_on);
01332 qc_lampflat.std_of=cpl_vector_get_stdev(vec_of);
01333 qc_lampflat.std_di=cpl_vector_get_stdev(vec_di);
01334 }
01335 qc_lampflat.nsat=cpl_vector_get_mean(vec_nsat);
01336 cpl_vector_delete(vec_on);
01337 cpl_vector_delete(vec_of);
01338 cpl_vector_delete(vec_di);
01339 cpl_vector_delete(vec_nsat);
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349 return 0;
01350 }
01351