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 #ifdef HAVE_CONFIG_H
00031 #include <config.h>
00032 #endif
00033
00034 #include <cpl.h>
00035 #include <stdio.h>
00036 #include <math.h>
00037
00038 #include "vircam_utils.h"
00039 #include "vircam_mask.h"
00040 #include "vircam_pfits.h"
00041 #include "vircam_dfs.h"
00042 #include "vircam_stats.h"
00043 #include "vircam_paf.h"
00044
00045
00046
00047 static int vircam_dark_current_create(cpl_plugin *) ;
00048 static int vircam_dark_current_exec(cpl_plugin *) ;
00049 static int vircam_dark_current_destroy(cpl_plugin *) ;
00050 static int vircam_dark_current(cpl_parameterlist *, cpl_frameset *) ;
00051 static int vircam_dark_current_save(cpl_frameset *filelist,
00052 cpl_parameterlist *parlist);
00053 static int vircam_dark_current_lastbit(int jext, cpl_frameset *framelist,
00054 cpl_parameterlist *parlist);
00055 static void vircam_dark_current_init(void);
00056 static void vircam_dark_current_tidy();
00057
00058
00059
00060 static struct {
00061
00062
00063
00064 float thresh;
00065 int extenum;
00066
00067
00068
00069 float mean_dark_current;
00070 } vircam_dark_current_config;
00071
00072 static struct {
00073 int *labels;
00074 cpl_frameset *darklist;
00075 vir_mask *master_mask;
00076 int nframes;
00077 float **data;
00078 vir_fits **allfits;
00079 double *subset;
00080 double *exps;
00081 cpl_image *outimage;
00082 vir_fits **good;
00083 int ngood;
00084 cpl_propertylist *phupaf;
00085 } ps;
00086
00087 static cpl_frame *product_frame = NULL;
00088 static int isfirst;
00089 static int dummy;
00090
00091 static char vircam_dark_current_description[] =
00092 "vircam_dark_current -- VIRCAM recipe for measuring dark current.\n"
00093 "A list of dark frames is given. A robust estimate of the dark current\n"
00094 "is calculated by fitting a median slope to the dark value vs exposure time\n"
00095 "for each pixel. The output is to a dark current map which shows the dark\n"
00096 "current in counts per second for each input pixel.\n\n"
00097 "The program requires the following files in the SOF:\n\n"
00098 " Tag Description\n"
00099 " -----------------------------------------------------------------------\n"
00100 " %-21s A list of raw dark images with various exposure times\n"
00101 " %-21s Optional master bad pixel map or\n"
00102 " %-21s Optional master confidence map\n"
00103 "\n";
00104
00167
00168
00169
00177
00178
00179 int cpl_plugin_get_info(cpl_pluginlist *list) {
00180 cpl_recipe *recipe = cpl_calloc(1,sizeof(*recipe));
00181 cpl_plugin *plugin = &recipe->interface;
00182 char alldesc[SZ_ALLDESC];
00183 (void)snprintf(alldesc,SZ_ALLDESC,vircam_dark_current_description,
00184 VIRCAM_DARKCUR_RAW,VIRCAM_CAL_BPM,VIRCAM_CAL_CONF);
00185
00186 cpl_plugin_init(plugin,
00187 CPL_PLUGIN_API,
00188 VIRCAM_BINARY_VERSION,
00189 CPL_PLUGIN_TYPE_RECIPE,
00190 "vircam_dark_current",
00191 "VIRCAM recipe to determine detector dark current",
00192 alldesc,
00193 "Jim Lewis",
00194 "jrl@ast.cam.ac.uk",
00195 vircam_get_license(),
00196 vircam_dark_current_create,
00197 vircam_dark_current_exec,
00198 vircam_dark_current_destroy);
00199
00200 cpl_pluginlist_append(list,plugin);
00201
00202 return(0);
00203 }
00204
00205
00214
00215
00216 static int vircam_dark_current_create(cpl_plugin *plugin) {
00217 cpl_recipe *recipe;
00218 cpl_parameter *p;
00219
00220
00221
00222 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
00223 recipe = (cpl_recipe *)plugin;
00224 else
00225 return(-1);
00226
00227
00228
00229 recipe->parameters = cpl_parameterlist_new();
00230
00231
00232
00233 p = cpl_parameter_new_value("vircam.vircam_dark_current.thresh",
00234 CPL_TYPE_DOUBLE,
00235 "Rejection threshold in sigma above background", "vircam.vircam_dark_current",5.0);
00236 cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"thresh");
00237 cpl_parameterlist_append(recipe->parameters,p);
00238
00239
00240
00241 p = cpl_parameter_new_range("vircam.vircam_dark_current.extenum",
00242 CPL_TYPE_INT,
00243 "Extension number to be done, 0 == all",
00244 "vircam.vircam_dark_current",
00245 1,0,16);
00246 cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"ext");
00247 cpl_parameterlist_append(recipe->parameters,p);
00248
00249
00250
00251 return(0);
00252 }
00253
00254
00255
00261
00262
00263 static int vircam_dark_current_destroy(cpl_plugin *plugin) {
00264 cpl_recipe *recipe ;
00265
00266
00267
00268 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
00269 recipe = (cpl_recipe *)plugin;
00270 else
00271 return(-1);
00272
00273 cpl_parameterlist_delete(recipe->parameters);
00274 return(0);
00275 }
00276
00277
00278
00284
00285
00286 static int vircam_dark_current_exec(cpl_plugin *plugin) {
00287 cpl_recipe *recipe;
00288
00289
00290
00291 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
00292 recipe = (cpl_recipe *)plugin;
00293 else
00294 return(-1);
00295
00296 return(vircam_dark_current(recipe->parameters,recipe->frames));
00297 }
00298
00299
00306
00307
00308 static int vircam_dark_current(cpl_parameterlist *parlist,
00309 cpl_frameset *framelist) {
00310 int jst,jfn,nlab,i,j,nx,ny,n,retval,live;
00311 long npts;
00312 double intercept,slope,sig;
00313 const char *fctid = "vircam_dark_current";
00314 float *outdata,val;
00315 unsigned char *bpm;
00316 vir_fits *ff;
00317 cpl_frame *cur_frame;
00318 cpl_propertylist *plist;
00319 cpl_parameter *p;
00320
00321
00322
00323 if (framelist == NULL || cpl_frameset_get_size(framelist) <= 0) {
00324 cpl_msg_error(fctid,"Input framelist NULL or has no input data\n");
00325 return(-1);
00326 }
00327
00328
00329
00330 vircam_dark_current_init();
00331
00332
00333
00334 p = cpl_parameterlist_find(parlist,"vircam.vircam_dark_current.thresh");
00335 vircam_dark_current_config.thresh = (float)cpl_parameter_get_double(p);
00336 p = cpl_parameterlist_find(parlist,"vircam.vircam_dark_current.extenum");
00337 vircam_dark_current_config.extenum = cpl_parameter_get_int(p);
00338
00339
00340
00341 if (vircam_dfs_set_groups(framelist) != VIR_OK) {
00342 cpl_msg_error(fctid,"Cannot identify RAW and CALIB frames");
00343 return(-1);
00344 }
00345
00346
00347
00348 if ((ps.labels = cpl_frameset_labelise(framelist,vircam_compare_tags,
00349 &nlab)) == NULL) {
00350 cpl_msg_error(fctid,"Cannot labelise the input frameset");
00351 vircam_dark_current_tidy();
00352 return(-1);
00353 }
00354 if ((ps.darklist = vircam_frameset_subgroup(framelist,ps.labels,nlab,
00355 VIRCAM_DARKCUR_RAW)) == NULL) {
00356 cpl_msg_error(fctid,"Cannot find dark frames in input frameset");
00357 vircam_dark_current_tidy();
00358 return(-1);
00359 }
00360 ps.nframes = cpl_frameset_get_size(ps.darklist);
00361 if (ps.nframes < 2) {
00362 cpl_msg_error(fctid,"Dark frameset doesn't have enough frames");
00363 vircam_dark_current_tidy();
00364 return(-1);
00365 }
00366
00367
00368
00369
00370 ps.master_mask = vircam_mask_define(framelist,ps.labels,nlab);
00371
00372
00373
00374 ps.data = cpl_malloc(ps.nframes*sizeof(float *));
00375 ps.subset = cpl_malloc(ps.nframes*sizeof(double));
00376 ps.exps = cpl_malloc(ps.nframes*sizeof(double));
00377
00378
00379
00380 for (i = 0; i < ps.nframes; i++) {
00381 cur_frame = cpl_frameset_get_frame(ps.darklist,i);
00382 plist = cpl_propertylist_load(cpl_frame_get_filename(cur_frame),0);
00383 if (vircam_pfits_get_exptime(plist,&val) != VIR_OK) {
00384 cpl_msg_error(fctid,"Unable to get exposure time for %s",
00385 cpl_frame_get_filename(cur_frame));
00386 return(-1);
00387 }
00388 ps.exps[i] = (double)val;
00389 cpl_propertylist_delete(plist);
00390 }
00391
00392
00393
00394
00395
00396 vircam_exten_range(vircam_dark_current_config.extenum,
00397 (const cpl_frame *)cpl_frameset_get_frame(ps.darklist,0),
00398 &jst,&jfn);
00399 if (jst == -1 || jfn == -1) {
00400 cpl_msg_error(fctid,"Unable to continue");
00401 vircam_dark_current_tidy();
00402 return(-1);
00403 }
00404
00405
00406
00407 ps.good = cpl_malloc(ps.nframes*sizeof(vir_fits *));
00408
00409
00410
00411 for (j = jst; j <= jfn; j++) {
00412 isfirst = (j == jst);
00413 dummy = 0;
00414 vircam_dark_current_config.mean_dark_current = 0.0;
00415
00416
00417
00418
00419 ps.allfits = vircam_fits_load_list(ps.darklist,CPL_TYPE_FLOAT,j);
00420 if (ps.allfits == NULL) {
00421 cpl_msg_info(fctid,"Extension %d darks wouldn't load",j);
00422 dummy = 1;
00423 retval = vircam_dark_current_lastbit(j,framelist,parlist);
00424 if (retval != 0)
00425 return(-1);
00426 continue;
00427 }
00428
00429
00430
00431 ps.ngood = 0;
00432 for (i = 0; i < ps.nframes; i++) {
00433 ff = ps.allfits[i];
00434 vircam_pfits_get_detlive(vircam_fits_get_ehu(ff),&live);
00435 if (! live) {
00436 cpl_msg_info(fctid,"Detector flagged dead %s",
00437 vircam_fits_get_fullname(ff));
00438 vircam_fits_set_error(ff,VIR_FATAL);
00439 } else {
00440 ps.good[ps.ngood] = ff;
00441 ps.ngood += 1;
00442 }
00443 }
00444
00445
00446
00447
00448 if (ps.ngood < 2) {
00449 cpl_msg_warning(fctid,"Need at least 2 good images -- %d found",
00450 ps.ngood);
00451 dummy = 1;
00452 retval = vircam_dark_current_lastbit(j,framelist,parlist);
00453 freefitslist(ps.allfits,ps.nframes);
00454 freeimage(ps.outimage);
00455 if (retval != 0)
00456 return(-1);
00457 continue;
00458 }
00459
00460
00461
00462 for (i = 0; i < ps.ngood; i++)
00463 ps.data[i] = cpl_image_get_data(vircam_fits_get_image(ps.allfits[i]));
00464
00465
00466
00467 nx = cpl_image_get_size_x(vircam_fits_get_image(ps.good[0]));
00468 ny = cpl_image_get_size_y(vircam_fits_get_image(ps.good[0]));
00469 retval = vircam_mask_load(ps.master_mask,j,nx,ny);
00470 if (retval == VIR_FATAL) {
00471 cpl_msg_info(fctid,"Unable to load mask image %s[%d]",
00472 vircam_mask_get_filename(ps.master_mask),j);
00473 cpl_msg_info(fctid,"Forcing all pixels to be good from now on");
00474 vircam_mask_force(ps.master_mask,nx,ny);
00475 }
00476 bpm = vircam_mask_get_data(ps.master_mask);
00477
00478
00479
00480 ps.outimage = cpl_image_new(nx,ny,CPL_TYPE_FLOAT);
00481 outdata = cpl_image_get_data(ps.outimage);
00482
00483
00484
00485 cpl_msg_info(fctid,"Doing dark current fits for extension %d",j);
00486 npts = (long)(nx*ny);
00487 for (n = 0; n < npts; n++) {
00488 if (bpm[n] != 0) {
00489 slope = 0.0;
00490 } else {
00491 for (i = 0; i < ps.ngood; i++)
00492 ps.subset[i] = (double)(ps.data[i][n]);
00493 vircam_linfit(ps.ngood,ps.exps,ps.subset,&intercept,&slope,
00494 &sig);
00495 }
00496
00497
00498
00499 outdata[n] = (float)slope;
00500 }
00501
00502
00503
00504 vircam_dark_current_config.mean_dark_current =
00505 vircam_med(outdata,bpm,npts);
00506
00507
00508
00509 (void)vircam_dark_current_lastbit(j,framelist,parlist);
00510
00511
00512
00513 freeimage(ps.outimage);
00514 vircam_mask_clear(ps.master_mask);
00515 freefitslist(ps.allfits,ps.nframes);
00516 }
00517
00518
00519
00520 vircam_dark_current_tidy();
00521
00522 return(0);
00523 }
00524
00525
00532
00533
00534 static int vircam_dark_current_save(cpl_frameset *framelist,
00535 cpl_parameterlist *parlist) {
00536 cpl_propertylist *plist,*p;
00537 const char *fctid = "vircam_dark_current_save";
00538 const char *outfile = "darkcurrent.fits";
00539 const char *outpaf = "darkcurrent";
00540 const char *recipeid = "vircam_dark_current";
00541 float darkcur_med;
00542
00543
00544
00545
00546 darkcur_med = vircam_dark_current_config.mean_dark_current;
00547 if (isfirst) {
00548
00549
00550
00551 product_frame = cpl_frame_new();
00552 cpl_frame_set_filename(product_frame,outfile);
00553 cpl_frame_set_tag(product_frame,VIRCAM_PRO_DARKCUR);
00554 cpl_frame_set_type(product_frame,CPL_FRAME_TYPE_IMAGE);
00555 cpl_frame_set_group(product_frame,CPL_FRAME_GROUP_PRODUCT);
00556 cpl_frame_set_level(product_frame,CPL_FRAME_LEVEL_FINAL);
00557
00558
00559
00560 plist = vircam_fits_get_phu(ps.allfits[0]);
00561 vircam_dfs_set_product_primary_header(plist,product_frame,framelist,
00562 parlist,(char *)recipeid,
00563 "PRO-1.15");
00564 ps.phupaf = vircam_paf_phu_items(plist);
00565
00566
00567
00568 if (cpl_image_save(NULL,outfile,CPL_BPP_8_UNSIGNED,plist,
00569 CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
00570 cpl_msg_error(fctid,"Cannot save product PHU");
00571 cpl_propertylist_delete(plist);
00572 return(-1);
00573 }
00574 cpl_frameset_insert(framelist,product_frame);
00575 }
00576
00577
00578
00579 plist = vircam_fits_get_ehu(ps.allfits[0]);
00580 vircam_dfs_set_product_exten_header(plist,product_frame,framelist,
00581 parlist,(char *)recipeid,
00582 "PRO-1.15");
00583
00584
00585
00586 cpl_propertylist_update_float(plist,"ESO QC DARKCURRENT",darkcur_med);
00587 cpl_propertylist_set_comment(plist,"ESO QC DARKCURRENT",
00588 "[ADU/s] Median dark current");
00589 if (dummy)
00590 vircam_dummy_property(plist);
00591
00592
00593
00594 if (cpl_image_save(ps.outimage,outfile,CPL_BPP_IEEE_FLOAT,plist,
00595 CPL_IO_EXTEND) != CPL_ERROR_NONE) {
00596 cpl_msg_error(fctid,"Cannot save product image extension");
00597 cpl_propertylist_delete(plist);
00598 return(-1);
00599 }
00600
00601
00602
00603 p = vircam_paf_req_items(plist);
00604 vircam_merge_propertylists(p,ps.phupaf);
00605 if (vircam_paf_print((char *)outpaf,"VIRCAM/vircam_dark_current","QC file",
00606 p) != VIR_OK)
00607 cpl_msg_warning(fctid,"Unable to write PAF\n");
00608 cpl_propertylist_delete(p);
00609
00610
00611
00612 return(0);
00613 }
00614
00615
00623
00624
00625 static int vircam_dark_current_lastbit(int jext, cpl_frameset *framelist,
00626 cpl_parameterlist *parlist) {
00627 int retval;
00628 const char *fctid="vircam_dark_current_lastbit";
00629
00630
00631
00632 if (dummy)
00633 ps.outimage = vircam_dummy_image(ps.allfits[0]);
00634
00635
00636
00637 cpl_msg_info(fctid,"Saving products for extension %d",jext);
00638 retval = vircam_dark_current_save(framelist,parlist);
00639 if (retval != 0) {
00640 vircam_dark_current_tidy();
00641 return(-1);
00642 }
00643 return(0);
00644 }
00645
00646
00647
00651
00652
00653 static void vircam_dark_current_init(void) {
00654 ps.labels = NULL;
00655 ps.darklist = NULL;
00656 ps.master_mask = NULL;
00657 ps.data = NULL;
00658 ps.allfits = NULL;
00659 ps.subset = NULL;
00660 ps.exps = NULL;
00661 ps.outimage = NULL;
00662 ps.nframes = 0;
00663 ps.good = NULL;
00664 ps.phupaf = NULL;
00665 }
00666
00667
00671
00672
00673 static void vircam_dark_current_tidy(void) {
00674
00675 freespace(ps.labels);
00676 freeframeset(ps.darklist);
00677 freemask(ps.master_mask);
00678 freespace(ps.data);
00679 freespace(ps.subset);
00680 freespace(ps.exps);
00681 freeimage(ps.outimage);
00682 freefitslist(ps.allfits,ps.nframes);
00683 ps.nframes = 0;
00684 freespace(ps.good);
00685 freepropertylist(ps.phupaf);
00686 }
00687
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839