00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031
00032
00033
00034
00035
00036 #include <math.h>
00037 #include <cpl.h>
00038
00039 #include "irplib_calib.h"
00040
00041
00042 static int
00043 irplib_get_clean_mean_window(cpl_image* img,
00044 const int llx,
00045 const int lly,
00046 const int urx, int ury,
00047 const int kappa,
00048 const int nclip,
00049 double* clean_mean,
00050 double* clean_stdev);
00051
00052 static double irplib_pfits_get_dit(const cpl_propertylist * plist);
00053 static double irplib_pfits_get_exp_time(const cpl_propertylist* plist);
00054
00058
00059
00061
00068
00069 static double irplib_pfits_get_dit(const cpl_propertylist * plist)
00070 {
00071 return cpl_propertylist_get_double(plist,"ESO DET DIT");
00072 }
00073
00074
00080
00081 static double irplib_pfits_get_exp_time(const cpl_propertylist* plist)
00082 {
00083
00084 return cpl_propertylist_get_double(plist,"EXPTIME");
00085
00086 }
00087
00088
00104 static int
00105 irplib_get_clean_mean_window(cpl_image* img,
00106 const int llx,
00107 const int lly,
00108 const int urx, int ury,
00109 const int kappa,
00110 const int nclip,
00111 double* clean_mean,
00112 double* clean_stdev)
00113 {
00114
00115
00116 double mean=0;
00117 double stdev=0;
00118 double threshold=0;
00119 double lo_cut=0;
00120 double hi_cut=0;
00121 cpl_mask* mask=NULL;
00122 cpl_image* tmp=NULL;
00123 cpl_stats* stats=NULL;
00124 int i=0;
00125
00126 tmp=cpl_image_extract(img,llx,lly,urx,ury);
00127 cpl_image_accept_all(tmp);
00128 for(i=0;i<nclip;i++) {
00129
00130
00131 cpl_stats_delete(stats);
00132 stats = cpl_stats_new_from_image(tmp, CPL_STATS_MEAN | CPL_STATS_STDEV);
00133 mean = cpl_stats_get_mean(stats);
00134 stdev = cpl_stats_get_stdev(stats);
00135
00136 threshold=kappa*stdev;
00137 lo_cut=mean-threshold;
00138 hi_cut=mean+threshold;
00139
00140 cpl_image_accept_all(tmp);
00141 mask=cpl_mask_threshold_image_create(tmp,lo_cut,hi_cut);
00142
00143 cpl_mask_not(mask);
00144 cpl_image_reject_from_mask(tmp,mask);
00145 cpl_mask_delete(mask);
00146
00147
00148 }
00149 *clean_mean=mean;
00150 *clean_stdev=stdev;
00151 cpl_image_delete(tmp);
00152 cpl_stats_delete(stats);
00153
00154 return 0;
00155
00156
00157 }
00158
00159
00176
00177
00178
00179 cpl_table*
00180 irplib_compute_gain(
00181 cpl_frameset* son,
00182 cpl_frameset* sof,
00183 int* zone1,
00184 int* zone2,
00185 const int kappa,
00186 const int nclip)
00187 {
00188
00189 cpl_frame* frm=NULL;
00190
00191 cpl_image* img_on1=NULL;
00192 cpl_image* img_on2=NULL;
00193 cpl_image* img_on_dif=NULL;
00194
00195 cpl_image* img_of1=NULL;
00196 cpl_image* img_of2=NULL;
00197 cpl_image* img_of_dif=NULL;
00198
00199 cpl_table* res_tbl=NULL;
00200 cpl_vector* dit_on=NULL;
00201 cpl_vector* dit_of=NULL;
00202 cpl_vector* exptime_on=NULL;
00203 cpl_vector* exptime_of=NULL;
00204 cpl_propertylist* plist=NULL;
00205
00206 int non=0;
00207 int nof=0;
00208 int nfr=0;
00209 int llx;
00210 int lly;
00211 int urx;
00212 int ury;
00213
00214 double avg_on1=0;
00215 double avg_on2=0;
00216 double avg_of1=0;
00217 double avg_of2=0;
00218 double std=0;
00219
00220 double sig_on_dif=0;
00221 double sig_of_dif=0;
00222 char* name=NULL;
00223 int i=0;
00224 int m=0;
00225
00226 double gain=0;
00227 double dit_ref=0;
00228 double dit_tmp=0;
00229 double exptime_ref=0;
00230 double exptime_tmp=0;
00231
00232
00233 non = cpl_frameset_get_size(son);
00234 nof = cpl_frameset_get_size(sof);
00235 nfr = (non <= nof) ? non : nof;
00236
00237 dit_on=cpl_vector_new(nfr);
00238 dit_of=cpl_vector_new(nfr);
00239 exptime_on=cpl_vector_new(nfr);
00240 exptime_of=cpl_vector_new(nfr);
00241
00242 for(i=0;i<nfr;i++) {
00243
00244 frm=cpl_frameset_get_frame(son,i);
00245 name=(char*)cpl_frame_get_filename(frm);
00246 plist=cpl_propertylist_load(name,0);
00247 dit_ref=irplib_pfits_get_dit(plist);
00248 exptime_ref=(double)irplib_pfits_get_exp_time(plist);
00249 cpl_propertylist_delete(plist);
00250 cpl_vector_set(dit_on,i,dit_ref);
00251 cpl_vector_set(exptime_on,i,exptime_ref);
00252
00253 frm=cpl_frameset_get_frame(sof,i);
00254 name=(char*)cpl_frame_get_filename(frm);
00255 plist=cpl_propertylist_load(name,0);
00256 dit_ref=irplib_pfits_get_dit(plist);
00257 exptime_ref=(double)irplib_pfits_get_exp_time(plist);
00258 cpl_propertylist_delete(plist);
00259 cpl_vector_set(dit_of,i,dit_ref);
00260 cpl_vector_set(exptime_of,i,exptime_ref);
00261
00262 }
00263
00264
00265 llx=zone1[0];
00266 lly=zone1[1];
00267 urx=zone1[2];
00268 ury=zone1[3];
00269
00270
00271
00272 res_tbl=cpl_table_new(nfr);
00273 cpl_table_new_column(res_tbl,"adu", CPL_TYPE_DOUBLE);
00274 cpl_table_new_column(res_tbl,"gain", CPL_TYPE_DOUBLE);
00275
00276 for(i=0;i<nfr;i++) {
00277 frm=cpl_frameset_get_frame(son,i);
00278 name=(char*)cpl_frame_get_filename(frm);
00279 img_on1=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00280
00281 frm=cpl_frameset_get_frame(sof,i);
00282 name=(char*)cpl_frame_get_filename(frm);
00283 img_of1=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00284
00285
00286 dit_ref=cpl_vector_get(dit_on,i);
00287 exptime_ref=cpl_vector_get(exptime_on,i);
00288
00289
00290 for(m=0;m<nfr; m++) {
00291 if(m != i) {
00292 frm=cpl_frameset_get_frame(son,m);
00293 name=(char*)cpl_frame_get_filename(frm);
00294 dit_tmp=cpl_vector_get(dit_on,m);
00295 exptime_tmp=cpl_vector_get(exptime_on,m);
00296 if(dit_tmp == dit_ref && exptime_tmp == exptime_ref) {
00297 img_on2=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00298 frm=cpl_frameset_get_frame(sof,m);
00299 name=(char*)cpl_frame_get_filename(frm);
00300 img_of2=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00301
00302 img_on_dif=cpl_image_subtract_create(img_on1,img_on2);
00303 img_of_dif=cpl_image_subtract_create(img_of1,img_of2);
00304
00305 irplib_get_clean_mean_window(img_on1,llx,lly,urx,ury,kappa,
00306 nclip,&avg_on1,&std);
00307 irplib_get_clean_mean_window(img_on2,llx,lly,urx,ury,kappa,
00308 nclip,&avg_on2,&std);
00309 irplib_get_clean_mean_window(img_of1,llx,lly,urx,ury,kappa,
00310 nclip,&avg_of1,&std);
00311 irplib_get_clean_mean_window(img_of2,llx,lly,urx,ury,kappa,
00312 nclip,&avg_of2,&std);
00313
00314 cpl_flux_get_noise_window(img_on_dif,zone2,2,100, &sig_on_dif, NULL);
00315 cpl_flux_get_noise_window(img_of_dif,zone2,2,100, &sig_of_dif, NULL);
00316
00317 cpl_image_delete(img_on2);
00318 cpl_image_delete(img_of2);
00319 cpl_image_delete(img_on_dif);
00320 cpl_image_delete(img_of_dif);
00321
00322 gain=((avg_on1+avg_on2)-(avg_of1+avg_of2))/
00323 ((sig_on_dif*sig_on_dif)-(sig_of_dif*sig_of_dif));
00324
00325 cpl_table_set_double(res_tbl,"gain",m,gain);
00326 cpl_table_set_double(res_tbl,"adu",m,
00327 ((avg_on1+avg_on2)/2-(avg_of1+avg_of2)/2));
00328
00329 }
00330 }
00331 }
00332 cpl_image_delete(img_on1);
00333 cpl_image_delete(img_of1);
00334 }
00335
00336
00337 cpl_vector_delete(dit_on);
00338 cpl_vector_delete(dit_of);
00339 cpl_vector_delete(exptime_on);
00340 cpl_vector_delete(exptime_of);
00341
00342 return res_tbl;
00343
00344 }
00345
00346
00356
00357
00358
00359 cpl_table* irplib_compute_linearity(cpl_frameset* son, cpl_frameset* sof)
00360 {
00361
00362
00363 cpl_frame* frm=NULL;
00364
00365 int* status=0;
00366 int non=0;
00367 int nof=0;
00368 int nfr=0;
00369 int i=0;
00370 double med_on=0;
00371 double avg_on=0;
00372 double med_of=0;
00373 double avg_of=0;
00374 double med_dit=0;
00375 double avg_dit=0;
00376
00377 double med=0;
00378 double avg=0;
00379
00380 char* name=NULL;
00381 cpl_image* img=NULL;
00382 cpl_vector* vec_adl=NULL;
00383 cpl_vector* vec_dit=NULL;
00384 cpl_vector* vec_avg=NULL;
00385 cpl_vector* vec_med=NULL;
00386 cpl_vector* vec_avg_dit=NULL;
00387 cpl_vector* vec_med_dit=NULL;
00388 cpl_propertylist* plist=NULL;
00389
00390 double dit=0;
00391 cpl_table* lin_tbl=NULL;
00392
00393
00394 non = cpl_frameset_get_size(son);
00395 nof = cpl_frameset_get_size(sof);
00396 nfr = (non <= nof) ? non : nof;
00397
00398 lin_tbl=cpl_table_new(nfr);
00399 cpl_table_new_column(lin_tbl,"med", CPL_TYPE_DOUBLE);
00400 cpl_table_new_column(lin_tbl,"avg", CPL_TYPE_DOUBLE);
00401 cpl_table_new_column(lin_tbl,"med_dit", CPL_TYPE_DOUBLE);
00402 cpl_table_new_column(lin_tbl,"avg_dit", CPL_TYPE_DOUBLE);
00403 cpl_table_new_column(lin_tbl,"dit", CPL_TYPE_DOUBLE);
00404 vec_med=cpl_vector_new(nfr);
00405 vec_avg=cpl_vector_new(nfr);
00406 vec_med_dit=cpl_vector_new(nfr);
00407 vec_avg_dit=cpl_vector_new(nfr);
00408 vec_dit=cpl_vector_new(nfr);
00409 vec_adl=cpl_vector_new(nfr);
00410 for(i=0;i<nfr;i++) {
00411 frm=cpl_frameset_get_frame(son,i);
00412 name=(char*)cpl_frame_get_filename(frm);
00413 img=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00414 med_on=cpl_image_get_median(img);
00415 avg_on=cpl_image_get_mean(img);
00416 cpl_image_delete(img);
00417
00418 frm=cpl_frameset_get_frame(sof,i);
00419 name=(char*)cpl_frame_get_filename(frm);
00420 img=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00421 med_of=cpl_image_get_median(img);
00422 avg_of=cpl_image_get_mean(img);
00423 cpl_image_delete(img);
00424 med=med_on-med_of;
00425 avg=avg_on-avg_of;
00426 plist=cpl_propertylist_load(name,0);
00427 dit=(double)irplib_pfits_get_dit(plist);
00428 cpl_propertylist_delete(plist);
00429 avg_dit=avg/dit;
00430 med_dit=med/dit;
00431
00432 cpl_vector_set(vec_dit,i,dit);
00433 cpl_vector_set(vec_avg,i,avg);
00434 cpl_vector_set(vec_med,i,med);
00435 cpl_vector_set(vec_avg_dit,i,avg_dit);
00436 cpl_vector_set(vec_med_dit,i,med_dit);
00437
00438
00439 cpl_table_set_double(lin_tbl,"dit",i,dit);
00440 cpl_table_set_double(lin_tbl,"med",i,med);
00441 cpl_table_set_double(lin_tbl,"avg",i,avg);
00442 cpl_table_set_double(lin_tbl,"med_dit",i,med_dit);
00443 cpl_table_set_double(lin_tbl,"avg_dit",i,avg_dit);
00444
00445 }
00446 cpl_table_new_column(lin_tbl,"adl", CPL_TYPE_DOUBLE);
00447 med_dit=cpl_vector_get_mean(vec_med_dit);
00448 avg_dit=cpl_vector_get_mean(vec_avg_dit);
00449
00450 for(i=0;i<nfr;i++) {
00451 dit = cpl_table_get_double(lin_tbl,"dit",i,status);
00452 cpl_vector_set(vec_adl,i,dit*med_dit);
00453 cpl_table_set_double(lin_tbl,"adl",i,dit*med_dit);
00454 }
00455
00456
00457 cpl_vector_delete(vec_dit);
00458 cpl_vector_delete(vec_adl);
00459 cpl_vector_delete(vec_avg);
00460 cpl_vector_delete(vec_med);
00461 cpl_vector_delete(vec_avg_dit);
00462 cpl_vector_delete(vec_med_dit);
00463
00464
00465 return lin_tbl;
00466
00467 }
00468
00469
00470
00479
00480 int irplib_detlin_correct(
00481 cpl_imagelist * ilist,
00482 const char * detlin_a,
00483 const char * detlin_b,
00484 const char * detlin_c)
00485 {
00486 cpl_image * ima ;
00487 cpl_image * imb ;
00488 cpl_image * imc ;
00489 float * pima ;
00490 float * pimb ;
00491 float * pimc ;
00492 float * pdata ;
00493 int nx, ny, ni ;
00494 double coeff_1, coeff_2, val ;
00495 int i, j ;
00496
00497
00498 if (!ilist || !detlin_a || !detlin_b || !detlin_c) return -1 ;
00499
00500
00501 ima = cpl_image_load(detlin_a, CPL_TYPE_FLOAT, 0, 0) ;
00502 imb = cpl_image_load(detlin_b, CPL_TYPE_FLOAT, 0, 0) ;
00503 imc = cpl_image_load(detlin_c, CPL_TYPE_FLOAT, 0, 0) ;
00504 if (!ima || !imb || !imc) {
00505 cpl_msg_error(cpl_func, "Cannot load the detlin images") ;
00506 if (ima) cpl_image_delete(ima) ;
00507 if (imb) cpl_image_delete(imb) ;
00508 if (imc) cpl_image_delete(imc) ;
00509 return -1 ;
00510 }
00511 pima = cpl_image_get_data_float(ima) ;
00512 pimb = cpl_image_get_data_float(imb) ;
00513 pimc = cpl_image_get_data_float(imc) ;
00514
00515
00516 nx = cpl_image_get_size_x(cpl_imagelist_get(ilist, 0)) ;
00517 ny = cpl_image_get_size_y(cpl_imagelist_get(ilist, 0)) ;
00518 ni = cpl_imagelist_get_size(ilist) ;
00519 if ((cpl_image_get_size_x(ima) != nx) ||
00520 (cpl_image_get_size_x(imb) != nx) ||
00521 (cpl_image_get_size_x(imc) != nx) ||
00522 (cpl_image_get_size_y(ima) != ny) ||
00523 (cpl_image_get_size_y(imb) != ny) ||
00524 (cpl_image_get_size_y(imc) != ny)) {
00525 cpl_msg_error(cpl_func, "Incompatible sizes") ;
00526 cpl_image_delete(ima) ;
00527 cpl_image_delete(imb) ;
00528 cpl_image_delete(imc) ;
00529 return -1 ;
00530 }
00531
00532
00533 for (i=0 ; i<nx*ny ; i++) {
00534
00535 if (fabs(pima[i]) < 1e-30) {
00536 coeff_1 = coeff_2 = (double)0.0 ;
00537 } else {
00538 coeff_1 = (double)pimb[i] / (double)pima[i] ;
00539 coeff_2 = (double)pimc[i] / (double)pima[i] ;
00540 }
00541
00542 for (j=0 ; j<ni ; j++) {
00543 pdata = cpl_image_get_data_float(cpl_imagelist_get(ilist, j)) ;
00544 val = (double)pdata[i] ;
00545 pdata[i]=(float)(val+coeff_1*val*val+coeff_2*val*val*val) ;
00546 }
00547 }
00548
00549 cpl_image_delete(ima) ;
00550 cpl_image_delete(imb) ;
00551 cpl_image_delete(imc) ;
00552 return 0 ;
00553 }
00554
00555
00564
00565 int irplib_flat_dark_bpm_calib(
00566 cpl_imagelist * ilist,
00567 const char * flat,
00568 const char * dark,
00569 const char * bpm)
00570 {
00571 cpl_image * dark_image ;
00572 cpl_image * flat_image ;
00573 cpl_mask * bpm_im_bin ;
00574 cpl_image * bpm_im_int ;
00575 int i ;
00576
00577
00578 if (ilist == NULL) return -1 ;
00579
00580
00581 if (dark != NULL) {
00582 cpl_msg_info(cpl_func, "Subtract the dark to the images") ;
00583
00584 if ((dark_image = cpl_image_load(dark, CPL_TYPE_FLOAT, 0, 0)) == NULL) {
00585 cpl_msg_error(cpl_func, "Cannot load the dark %s", dark) ;
00586 return -1 ;
00587 }
00588
00589 if (cpl_imagelist_subtract_image(ilist, dark_image)!=CPL_ERROR_NONE) {
00590 cpl_msg_error(cpl_func, "Cannot apply the dark to the images") ;
00591 cpl_image_delete(dark_image) ;
00592 return -1 ;
00593 }
00594 cpl_image_delete(dark_image) ;
00595 }
00596
00597
00598 if (flat != NULL) {
00599 cpl_msg_info(cpl_func, "Divide the images by the flatfield") ;
00600
00601 if ((flat_image = cpl_image_load(flat, CPL_TYPE_FLOAT, 0, 0)) == NULL) {
00602 cpl_msg_error(cpl_func, "Cannot load the flat field %s", flat) ;
00603 return -1 ;
00604 }
00605
00606 if (cpl_imagelist_divide_image(ilist, flat_image)!=CPL_ERROR_NONE) {
00607 cpl_msg_error(cpl_func, "Cannot apply the flatfield to the images") ;
00608 cpl_image_delete(flat_image) ;
00609 return -1 ;
00610 }
00611 cpl_image_delete(flat_image) ;
00612 }
00613
00614
00615 if (bpm != NULL) {
00616 cpl_msg_info(cpl_func, "Correct the bad pixels in the images") ;
00617
00618 if ((bpm_im_int = cpl_image_load(bpm, CPL_TYPE_INT, 0, 0)) == NULL) {
00619 cpl_msg_error(cpl_func, "Cannot load the bad pixel map %s", bpm) ;
00620 return -1 ;
00621 }
00622
00623 bpm_im_bin = cpl_mask_threshold_image_create(bpm_im_int, -0.5, 0.5) ;
00624 cpl_mask_not(bpm_im_bin) ;
00625 cpl_image_delete(bpm_im_int) ;
00626
00627 for (i=0 ; i<cpl_imagelist_get_size(ilist) ; i++) {
00628 cpl_image_reject_from_mask(cpl_imagelist_get(ilist, i), bpm_im_bin);
00629 if (cpl_detector_interpolate_rejected(
00630 cpl_imagelist_get(ilist, i)) != CPL_ERROR_NONE) {
00631 cpl_msg_error(cpl_func, "Cannot clean the bad pixels in obj %d",
00632 i+1);
00633 cpl_mask_delete(bpm_im_bin) ;
00634 return -1 ;
00635 }
00636 }
00637 cpl_mask_delete(bpm_im_bin) ;
00638 }
00639
00640
00641 return 0 ;
00642 }
00643