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
00032 #ifdef HAVE_CONFIG_H
00033 # include <config.h>
00034 #endif
00035
00036
00037
00038
00039 #include <irplib_utils.h>
00040 #include <irplib_error.h>
00041 #include <assert.h>
00042 #include "sinfo_new_psf.h"
00043 #include "sinfo_pro_save.h"
00044 #include "sinfo_key_names.h"
00045 #include "sinfo_psf_ini.h"
00046 #include "sinfo_psf_ini_by_cpl.h"
00047 #include "sinfo_utilities_scired.h"
00048 #include "sinfo_hidden.h"
00049 #include "sinfo_pfits.h"
00050 #include "sinfo_functions.h"
00051 #include "sinfo_error.h"
00052 #include "sinfo_utils_wrappers.h"
00053 #include "sinfo_globals.h"
00054
00055
00056
00057
00058
00059 #define SINFO_STREHL_M1 8.0
00060 #define SINFO_STREHL_M2 1.1
00061 #define SINFO_STREHL_BOX_SIZE 64
00062 #define SINFO_STREHL_WINDOW 6
00063 #define SINFO_PSF_SZ 3
00064 #define SINFO_PI 3.1415926535897932384626433
00065 #define SINFO_STREHL_ERROR_COEFFICIENT SINFO_PI * 0.007 / 0.0271
00066
00067
00068
00069 static void
00070 sinfo_check_borders(int* val,const int max,const int thresh);
00071
00072 static void
00073 sinfo_get_safe_box(int* llx,
00074 int* lly,
00075 int* urx,
00076 int* ury,
00077 const int xpos,
00078 const int ypos,
00079 const int box,
00080 const int szx,
00081 const int szy);
00082
00083 static int
00084 sinfo_get_strehl_from_slice(cpl_imagelist* cube,
00085 double disp,
00086 double cWave,
00087 double ws,
00088 double we,
00089 double pscale,
00090 double strehl_star_radius,
00091 double strehl_bg_r1,
00092 double strehl_bg_r2,
00093 double* strehl,
00094 double* strehl_err);
00095
00096
00097 static cpl_table*
00098 sinfo_get_encircled_energy(cpl_frameset* sof,
00099 cpl_image* img,
00100 double* fwhm_x,
00101 double* fwhm_y,
00102 cpl_table** qclog);
00103
00104 static double
00105 sinfo_get_strehl_from_ima(cpl_image* ima,
00106 char* name,
00107 cpl_frame* frame);
00108
00109 static int
00110 sinfo_get_strehl_from_image(cpl_image* img,
00111 double ws,
00112 double we,
00113 double pscale,
00114 double strehl_star_radius,
00115 double strehl_bg_r1,
00116 double strehl_bg_r2,
00117 double* strehl,
00118 double* strehl_err);
00119
00120 static cpl_table*
00121 sinfo_get_strehl_from_cube(cpl_imagelist* cube,
00122 char* name,
00123 cpl_frame* frame);
00124
00125 static int
00126 sinfo_get_frm12(cpl_frameset* sof,cpl_frame** frm1,cpl_frame** frm2);
00127
00128 static cpl_table*
00129 sinfo_get_strehl_from_2cubes(cpl_imagelist* cube1,
00130 cpl_imagelist* cube2,
00131 const char* fname1,
00132 cpl_frame* frm1,
00133 cpl_frame* frm2);
00134
00135 static int sinfo_strehl_compute(
00136 const cpl_image * im1,
00137 const cpl_image * im2,
00138 double m1,
00139 double m2,
00140 double lam,
00141 double dlam,
00142 double pscale1,
00143 double pscale2,
00144 int size,
00145 int xpos1,
00146 int ypos1,
00147 int xpos2,
00148 int ypos2,
00149 double r1,
00150 double r2,
00151 double r3,
00152 int noise_box_sz,
00153 int noise_nsamples,
00154 double * strehl,
00155 double * strehl_err,
00156 double * star_bg,
00157 double * star_peak,
00158 double * star_flux,
00159 double * psf_peak,
00160 double * psf_flux,
00161 double * bg_noise);
00162
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185 int
00186 sinfo_new_psf (const char* plugin_id,
00187 cpl_parameterlist* config,
00188 cpl_frameset* sof)
00189 {
00190
00191 cpl_imagelist* cube1=NULL;
00192 cpl_imagelist* cube2=NULL;
00193 cpl_image * med_img1=NULL ;
00194 cpl_image * med_img2=NULL ;
00195
00196 cpl_table* ao_performance=NULL;
00197 cpl_table* enc_energy=NULL;
00198
00199 cpl_frame* frm1=NULL;
00200 cpl_frame* frm2=NULL;
00201
00202 cpl_table* qclog_tbl=NULL;
00203 cpl_frameset* stk=NULL;
00204 cpl_propertylist* plist =NULL;
00205
00206 psf_config * cfg =NULL;
00207
00208 int nsample=0;
00209 int i = 0;
00210 int status=0;
00211
00212
00213
00214 int strehl_sw=0;
00215 int ilx1=0;
00216 int ily1=0;
00217 int ilx2=0;
00218 int ily2=0;
00219
00220 float cx1=0;
00221 float cy1=0;
00222 float cx2=0;
00223 float cy2=0;
00224
00225 double fwhm_x=0;
00226 double fwhm_y=0;
00227 double lam=0;
00228 double strehl=0;
00229
00230 char fname1[MAX_NAME_SIZE];
00231 char fname2[MAX_NAME_SIZE];
00232
00233 char key_name[MAX_NAME_SIZE];
00234
00235 char obs_name1[MAX_NAME_SIZE];
00236 char hlamp_st='F';
00237 char shut2_st='F';
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248 sinfo_msg("Parsing cpl input");
00249 check_nomsg(stk=cpl_frameset_new());
00250
00251 cknull(cfg = sinfo_parse_cpl_input_psf(sof,&stk),
00252 "error parsing cpl input");
00253
00254
00255 strehl_sw=sinfo_get_strehl_type(sof);
00256
00257 if(strehl_sw==0) {
00258 sinfo_msg("One target Strehl computation");
00259 if(sinfo_is_fits_file(cfg->inFrame) != 1) {
00260 sinfo_msg_error("Input file %s is not FITS",cfg->inFrame);
00261 goto cleanup;
00262 } else {
00263 strcpy(fname1,cfg->inFrame);
00264 }
00265
00266 if(NULL != cpl_frameset_find(sof,PRO_COADD_PSF)) {
00267 frm1 = cpl_frameset_find(sof,PRO_COADD_PSF);
00268 } else if(NULL != cpl_frameset_find(sof,PRO_OBS_PSF)) {
00269 frm1 = cpl_frameset_find(sof,PRO_OBS_PSF);
00270 } else if(NULL != cpl_frameset_find(sof,PRO_COADD_STD)) {
00271 frm1 = cpl_frameset_find(sof,PRO_COADD_STD);
00272 } else if(NULL != cpl_frameset_find(sof,PRO_OBS_STD)) {
00273 frm1 = cpl_frameset_find(sof,PRO_OBS_STD);
00274 } else if(NULL != cpl_frameset_find(sof,PRO_COADD_OBJ)) {
00275 frm1 = cpl_frameset_find(sof,PRO_COADD_OBJ);
00276 } else if(NULL != cpl_frameset_find(sof,PRO_OBS_OBJ)) {
00277 frm1 = cpl_frameset_find(sof,PRO_OBS_OBJ);
00278 } else {
00279 sinfo_msg_error("Frame %s or %s or %s or %s or %s or %s not found!",
00280 PRO_COADD_PSF,PRO_OBS_PSF,
00281 PRO_COADD_STD,PRO_OBS_STD,
00282 PRO_COADD_OBJ,PRO_OBS_OBJ);
00283 goto cleanup;
00284 }
00285
00286 sinfo_get_obsname(frm1,obs_name1);
00287 check_nomsg(hlamp_st=sinfo_get_keyvalue_bool(frm1,KEY_NAME_LAMP_HALO));
00288 check_nomsg(shut2_st=sinfo_get_keyvalue_bool(frm1,KEY_NAME_SHUT2_ST));
00289
00290 check_nomsg(cube1 = cpl_imagelist_load(fname1,CPL_TYPE_FLOAT,0));
00291 cknull(med_img1=sinfo_new_median_cube(cube1),
00292 " could not do sinfo_medianCube()");
00293
00294 check_nomsg(ilx1=cpl_image_get_size_x(med_img1));
00295 check_nomsg(ily1=cpl_image_get_size_y(med_img1));
00296
00297 cx1 = ilx1 / 2. + 0.5;
00298 cy1 = ily1 / 2. + 0.5;
00299 cknull(ao_performance=sinfo_get_strehl_from_cube(cube1,fname1,frm1),
00300 "error computing strehl");
00301 sinfo_free_imagelist(&cube1);
00302 } else {
00303 sinfo_msg("Two target Strehl computation");
00304 sinfo_get_frm12(sof,&frm1,&frm2);
00305 strcpy(fname1,cpl_frame_get_filename(frm1));
00306 strcpy(fname2,cpl_frame_get_filename(frm2));
00307
00308 check_nomsg(cube1 = cpl_imagelist_load(fname1,CPL_TYPE_FLOAT,0));
00309 check_nomsg(cube2 = cpl_imagelist_load(fname2,CPL_TYPE_FLOAT,0));
00310 cknull(med_img1=sinfo_new_median_cube(cube1),"Computing median on cube");
00311 cknull(med_img2=sinfo_new_median_cube(cube2),"Computing median on cube");
00312
00313
00314 check_nomsg(ilx1=cpl_image_get_size_x(med_img1));
00315 check_nomsg(ily1=cpl_image_get_size_y(med_img1));
00316 check_nomsg(ilx2=cpl_image_get_size_x(med_img2));
00317 check_nomsg(ily2=cpl_image_get_size_y(med_img2));
00318
00319 cx1 = ilx1 / 2. + 0.5;
00320 cy1 = ily1 / 2. + 0.5;
00321 cx2 = ilx2 / 2. + 0.5;
00322 cy2 = ily2 / 2. + 0.5;
00323 cknull(ao_performance=sinfo_get_strehl_from_2cubes(cube1,
00324 cube2,
00325 fname1,
00326 frm1,
00327 frm2),
00328 "Computing strehl");
00329
00330 sinfo_free_imagelist(&cube1);
00331 sinfo_free_imagelist(&cube2);
00332
00333 }
00334
00335
00336
00337 check_nomsg(nsample=cpl_table_get_nrow(ao_performance));
00338 cknull_nomsg(qclog_tbl = sinfo_qclog_init());
00339
00340
00341 cknull(plist = cpl_propertylist_load(fname1, 0),
00342 "getting header from reference ima frame %s",fname1);
00343
00344 if (cpl_propertylist_contains(plist, KEY_NAME_LOOP_STATE)) {
00345 sinfo_qclog_add_string(qclog_tbl,KEY_NAME_LOOP_STATE,
00346 cpl_propertylist_get_string(plist,KEY_NAME_LOOP_STATE),
00347 KEY_HELP_LOOP_STATE,"%s");
00348 }
00349
00350
00351 if (cpl_propertylist_contains(plist, KEY_NAME_LOOP_LGS)) {
00352 sinfo_qclog_add_int(qclog_tbl,KEY_NAME_LOOP_LGS,
00353 cpl_propertylist_get_int(plist,KEY_NAME_LOOP_LGS),
00354 KEY_HELP_LOOP_LGS,"%d");
00355 }
00356
00357 if (cpl_propertylist_contains(plist, KEY_NAME_INS1_MODE)) {
00358 sinfo_qclog_add_string(qclog_tbl,KEY_NAME_INS1_MODE,
00359 cpl_propertylist_get_string(plist,KEY_NAME_INS1_MODE),
00360 KEY_HELP_INS1_MODE,"%s");
00361 }
00362 sinfo_free_propertylist(&plist);
00363
00364 check_nomsg(strehl=cpl_table_get_column_median(ao_performance,"strehl"));
00365
00366 ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC STREHL MED",strehl,
00367 "STREHL MEDIAN","%f"));
00368 strehl=sinfo_get_strehl_from_ima(med_img1,fname1,frm1);
00369
00370 ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC STREHL AVG",strehl,
00371 "STREHL AVERAGE","%f"));
00372
00373 for(i=1;i<nsample;i++) {
00374
00375 check_nomsg(strehl=cpl_table_get_double(ao_performance,"strehl",
00376 i,&status));
00377 snprintf(key_name,MAX_NAME_SIZE-1,"%s%d","QC STREHL",i);
00378 ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,strehl,"STREHL","%f"));
00379
00380 check_nomsg(lam=cpl_table_get_double(ao_performance,"wavelength",
00381 i,&status));
00382 snprintf(key_name,MAX_NAME_SIZE-1,"%s%d","QC LAMBDA",i);
00383 ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,lam,
00384 "WAVELENGTH","%f"));
00385
00386 }
00387
00388 check_nomsg(strehl=cpl_table_get_column_median(ao_performance,
00389 "strehl_error"));
00390 ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC STREHL MEDERR",strehl,
00391 "STREHL ERROR MEDIAN","%f"));
00392 ck0_nomsg(sinfo_qclog_add_string(qclog_tbl,"OBS NAME",obs_name1,
00393 "OB name","%s"));
00394 ck0_nomsg(sinfo_qclog_add_bool(qclog_tbl,PAF_NAME_LAMP_HALO,hlamp_st,
00395 KEY_NAME_LAMP_HALO,"%d"));
00396 ck0_nomsg(sinfo_qclog_add_bool(qclog_tbl,PAF_NAME_SHUT2_ST,shut2_st,
00397 KEY_NAME_SHUT2_ST,"%d"));
00398
00399 ck0(sinfo_pro_save_tbl(ao_performance,stk,sof,
00400 PSF_AO_PERFORMANCE_OUT_FILENAME,
00401 PRO_AO_PERFORMANCE,qclog_tbl,plugin_id,config),
00402 "cannot save tbl %s", PSF_AO_PERFORMANCE_OUT_FILENAME);
00403
00404 sinfo_free_table(&qclog_tbl);
00405 sinfo_free_table(&ao_performance);
00406
00407
00408 cknull_nomsg(qclog_tbl=sinfo_qclog_init());
00409 cknull(enc_energy=sinfo_get_encircled_energy(sof,
00410 med_img1,
00411 &fwhm_x,
00412 &fwhm_y,
00413 &qclog_tbl),
00414 "Computing encircled energy");
00415
00416 ck0(sinfo_pro_save_tbl(enc_energy,stk,sof,PSF_ENC_ENERGY_OUT_FILENAME,
00417 PRO_ENC_ENERGY,qclog_tbl,plugin_id,config),
00418 "cannot save tbl %s", PSF_ENC_ENERGY_OUT_FILENAME);
00419
00420 sinfo_free_table(&qclog_tbl);
00421 sinfo_free_table(&enc_energy);
00422
00423
00424 cknull_nomsg(qclog_tbl = sinfo_qclog_init());
00425 ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC FWHMX",fwhm_x,
00426 "QC FWHM X","%f"));
00427 ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC FWHMY",fwhm_y,
00428 "QC FWHM Y","%f"));
00429 ck0_nomsg(sinfo_qclog_add_bool(qclog_tbl,PAF_NAME_LAMP_HALO,
00430 hlamp_st,KEY_NAME_LAMP_HALO,"%d"));
00431 ck0_nomsg(sinfo_qclog_add_bool(qclog_tbl,PAF_NAME_SHUT2_ST,shut2_st,
00432 KEY_NAME_SHUT2_ST,"%d"));
00433
00434 ck0(sinfo_pro_save_ima(med_img1,stk,sof,cfg->outName,PRO_PSF,
00435 qclog_tbl,plugin_id,config),
00436 "cannot save ima %s", cfg->outName);
00437
00438 sinfo_free_table(&qclog_tbl);
00439 sinfo_new_set_wcs_image(med_img1,cfg->outName,cx1, cy1);
00440 sinfo_free_image(&med_img1);
00441 sinfo_free_frameset(&stk);
00442 sinfo_free_psf(&cfg);
00443 return 0;
00444
00445 cleanup:
00446
00447 sinfo_free_table(&qclog_tbl);
00448 sinfo_free_imagelist(&cube2);
00449 sinfo_free_imagelist(&cube1);
00450 sinfo_free_table(&enc_energy);
00451 sinfo_free_image(&med_img1);
00452 sinfo_free_table(&ao_performance);
00453 sinfo_free_propertylist(&plist) ;
00454 sinfo_free_psf(&cfg);
00455 sinfo_free_frameset(&stk);
00456 return -1 ;
00457
00458 }
00459
00460 static int
00461 sinfo_get_strehl_from_image(cpl_image* img,
00462 double ws,
00463 double we,
00464 double pscale,
00465 double strehl_star_radius,
00466 double strehl_bg_r1,
00467 double strehl_bg_r2,
00468 double* strehl,
00469 double* strehl_err)
00470 {
00471
00472 cpl_image* img_dup=NULL;
00473
00474 double dlam=0.;
00475 double lam=0.;
00476 double norm=0.;
00477 double xc=0.;
00478 double yc=0.;
00479 double sx=0.;
00480 double sy=0.;
00481 double fwhm_x=0.;
00482 double fwhm_y=0.;
00483 double max_ima_cx=0.;
00484 double max_ima_cy=0.;
00485 double bkg=0.;
00486 double psf_peak=0.;
00487 double psf_flux=0.;
00488 double bkg_noise=0.;
00489 double star_bkg=0.;
00490 double star_peak=0.;
00491 double star_flux=0.;
00492
00493 int max_ima_x=0;
00494 int max_ima_y=0;
00495 int wllx=0;
00496 int wlly=0;
00497 int wurx=0;
00498 int wury=0;
00499 int ima_szx=0;
00500 int ima_szy=0;
00501
00502
00503 lam = (double)0.5*(ws+we);
00504 dlam=we-ws;
00505 sinfo_msg_warning("ws=%f we=%f dl=%f",ws,we,dlam);
00506 check_nomsg(img_dup=cpl_image_duplicate(img));
00507 sinfo_clean_nan(&img_dup);
00508 check_nomsg(cpl_image_get_maxpos(img_dup,&max_ima_x,&max_ima_y));
00509 sinfo_free_image(&img_dup);
00510
00511 check_nomsg(ima_szx=cpl_image_get_size_x(img));
00512 check_nomsg(ima_szy=cpl_image_get_size_y(img));
00513 sinfo_check_borders(&max_ima_x,ima_szx,SINFO_STREHL_WINDOW);
00514 sinfo_check_borders(&max_ima_y,ima_szy,SINFO_STREHL_WINDOW);
00515 sinfo_get_safe_box(&wllx,&wlly,&wurx,&wury,max_ima_x,max_ima_y,SINFO_PSF_SZ,
00516 ima_szx,ima_szy);
00517
00518
00519
00520
00521 check_nomsg(max_ima_cx=cpl_image_get_centroid_x_window(img,wllx,wlly,
00522 wurx,wury));
00523 check_nomsg(max_ima_cy=cpl_image_get_centroid_y_window(img,wllx,wlly,
00524 wurx,wury));
00525
00526 if(CPL_ERROR_NONE != cpl_image_fit_gaussian(img,max_ima_x,max_ima_y,
00527 SINFO_PSF_SZ,&norm,&xc,&yc,
00528 &sx,&sy,&fwhm_x,&fwhm_y)) {
00529
00530 sinfo_msg_warning("Gaussian fit failed");
00531 cpl_error_reset();
00532
00533 }
00534
00535
00536 check_nomsg(bkg=irplib_strehl_ring_background(img,max_ima_x,max_ima_y,
00537 19,21,IRPLIB_BG_METHOD_AVER_REJ)) ;
00538
00539
00540 if(CPL_ERROR_NONE != irplib_strehl_compute(img,
00541 SINFO_STREHL_M1,
00542 SINFO_STREHL_M2,
00543 lam,
00544 dlam,
00545 pscale,
00546 SINFO_STREHL_BOX_SIZE,
00547 max_ima_x,
00548 max_ima_y,
00549 strehl_star_radius,
00550 strehl_bg_r1,
00551 strehl_bg_r2,
00552 NOISE_HSIZE,
00553 NOISE_NSAMPLES,
00554 strehl,
00555 strehl_err,
00556 &star_bkg,
00557 &star_peak,
00558 &star_flux,
00559 &psf_peak,
00560 &psf_flux,
00561 &bkg_noise)) {
00562
00563
00564 sinfo_msg_warning("Problem computing strehl");
00565 *strehl=-1;
00566 *strehl_err=0;
00567
00568 irplib_error_reset();
00569 cpl_error_reset();
00570 }
00571
00572 return 0;
00573
00574 cleanup:
00575
00576 return -1;
00577
00578 }
00579
00580
00581
00582
00583 static int
00584 sinfo_get_strehl_from_slice(cpl_imagelist* cube,
00585 double disp,
00586 double cWave,
00587 double ws,
00588 double we,
00589 double pscale,
00590 double strehl_star_radius,
00591 double strehl_bg_r1,
00592 double strehl_bg_r2,
00593 double* strehl,
00594 double* strehl_err)
00595 {
00596 cpl_image* img_dup=NULL;
00597 cpl_image* img=NULL;
00598
00599 double dlam=0.;
00600 double lam=0.;
00601 double norm=0.;
00602 double xc=0.;
00603 double yc=0.;
00604 double sx=0.;
00605 double sy=0.;
00606 double fwhm_x=0.;
00607 double fwhm_y=0.;
00608 double max_ima_cx=0.;
00609 double max_ima_cy=0.;
00610 double bkg=0.;
00611 double psf_peak=0.;
00612 double psf_flux=0.;
00613 double bkg_noise=0.;
00614 double star_bkg=0.;
00615 double star_peak=0.;
00616 double star_flux=0.;
00617
00618 int max_ima_x=0;
00619 int max_ima_y=0;
00620 int wllx=0;
00621 int wlly=0;
00622 int wurx=0;
00623 int wury=0;
00624 int ima_szx=0;
00625 int ima_szy=0;
00626
00627
00628 lam = (double)0.5*(ws+we);
00629 dlam=we-ws;
00630
00631 img=sinfo_new_average_cube_to_image_between_waves(cube,disp,cWave,ws,we);
00632 check_nomsg(img_dup=cpl_image_duplicate(img));
00633 sinfo_clean_nan(&img_dup);
00634 check_nomsg(cpl_image_get_maxpos(img_dup,&max_ima_x,&max_ima_y));
00635 check_nomsg(cpl_image_delete(img_dup));
00636
00637 check_nomsg(ima_szx=cpl_image_get_size_x(img));
00638 check_nomsg(ima_szy=cpl_image_get_size_y(img));
00639 sinfo_check_borders(&max_ima_x,ima_szx,SINFO_STREHL_WINDOW);
00640 sinfo_check_borders(&max_ima_y,ima_szy,SINFO_STREHL_WINDOW);
00641 sinfo_get_safe_box(&wllx,&wlly,&wurx,&wury,max_ima_x,max_ima_y,SINFO_PSF_SZ,
00642 ima_szx,ima_szy);
00643
00644
00645
00646
00647 check_nomsg(max_ima_cx=cpl_image_get_centroid_x_window(img,wllx,wlly,
00648 wurx,wury));
00649 check_nomsg(max_ima_cy=cpl_image_get_centroid_y_window(img,wllx,wlly,
00650 wurx,wury));
00651
00652 if(CPL_ERROR_NONE != cpl_image_fit_gaussian(img,max_ima_x,max_ima_y,
00653 SINFO_PSF_SZ,&norm,&xc,&yc,
00654 &sx,&sy,&fwhm_x,&fwhm_y)) {
00655
00656 sinfo_msg_warning("Gaussian fit failed");
00657 cpl_error_reset();
00658
00659 }
00660
00661
00662 check_nomsg(bkg=irplib_strehl_ring_background(img,max_ima_x,max_ima_y,
00663 19,21,IRPLIB_BG_METHOD_AVER_REJ)) ;
00664
00665
00666 if(CPL_ERROR_NONE != irplib_strehl_compute(img,
00667 SINFO_STREHL_M1,
00668 SINFO_STREHL_M2,
00669 lam,
00670 dlam,
00671 pscale,
00672 SINFO_STREHL_BOX_SIZE,
00673 max_ima_x,
00674 max_ima_y,
00675 strehl_star_radius,
00676 strehl_bg_r1,
00677 strehl_bg_r2,
00678 NOISE_HSIZE,
00679 NOISE_NSAMPLES,
00680 strehl,
00681 strehl_err,
00682 &star_bkg,
00683 &star_peak,
00684 &star_flux,
00685 &psf_peak,
00686 &psf_flux,
00687 &bkg_noise)) {
00688
00689
00690 sinfo_msg_warning("Problem computing strehl");
00691 *strehl=-1;
00692 *strehl_err=0;
00693
00694 irplib_error_reset();
00695 cpl_error_reset();
00696 }
00697
00698
00699 sinfo_free_image(&img);
00700
00701
00702 return 0;
00703
00704 cleanup:
00705 return -1;
00706
00707 }
00708
00709
00710
00711
00712
00713 static int sinfo_get_strehl_from_2slices(cpl_imagelist* cube1,
00714 cpl_imagelist* cube2,
00715 double disp,
00716 double cWave,
00717 double ws,
00718 double we,
00719 double pscale1,
00720 double pscale2,
00721 double strehl_star_rad1,
00722 double strehl_bg_rmin1,
00723 double strehl_bg_rmax1,
00724 double* strehl,
00725 double* strehl_err)
00726 {
00727 cpl_image* img1=NULL;
00728 cpl_image* img2=NULL;
00729
00730 cpl_image* img_dup=NULL;
00731
00732 double dlam=0.;
00733 double lam=0.;
00734 double norm=0.;
00735 double xc=0.;
00736 double yc=0.;
00737 double sx=0.;
00738 double sy=0.;
00739 double fwhm_x=0.;
00740 double fwhm_y=0.;
00741 double max_ima1_cx=0.;
00742 double max_ima1_cy=0.;
00743
00744 double bkg1=0.;
00745 double psf_peak=0.;
00746 double psf_flux=0.;
00747 double bkg_noise=0.;
00748 double star_bkg=0.;
00749 double star_peak=0.;
00750 double star_flux=0.;
00751
00752 int max_ima1_x=0;
00753 int max_ima1_y=0;
00754
00755 int max_ima2_x=0;
00756 int max_ima2_y=0;
00757
00758 int wllx=0;
00759 int wlly=0;
00760 int wurx=0;
00761 int wury=0;
00762 int ima_szx=0;
00763 int ima_szy=0;
00764
00765
00766 lam = (double)0.5*(ws+we);
00767 dlam=we-ws;
00768
00769 img1=sinfo_new_average_cube_to_image_between_waves(cube1,disp,cWave,ws,we);
00770 img2=sinfo_new_average_cube_to_image_between_waves(cube2,disp,cWave,ws,we);
00771
00772 check_nomsg(img_dup=cpl_image_duplicate(img1));
00773 sinfo_clean_nan(&img_dup);
00774 check_nomsg(cpl_image_get_maxpos(img_dup,&max_ima1_x,&max_ima1_y));
00775 sinfo_free_image(&img_dup);
00776
00777 check_nomsg(img_dup=cpl_image_duplicate(img2));
00778 sinfo_clean_nan(&img_dup);
00779 check_nomsg(cpl_image_get_maxpos(img_dup,&max_ima2_x,&max_ima2_y));
00780 sinfo_free_image(&img_dup);
00781
00782 check_nomsg(ima_szx=cpl_image_get_size_x(img1));
00783 check_nomsg(ima_szy=cpl_image_get_size_y(img1));
00784 sinfo_check_borders(&max_ima1_x,ima_szx,SINFO_STREHL_WINDOW);
00785 sinfo_check_borders(&max_ima1_y,ima_szy,SINFO_STREHL_WINDOW);
00786 sinfo_get_safe_box(&wllx,&wlly,&wurx,&wury,max_ima1_x,max_ima1_y,
00787 SINFO_PSF_SZ,ima_szx,ima_szy);
00788
00789
00790
00791
00792 check_nomsg(max_ima1_cx=cpl_image_get_centroid_x_window(img1,wllx,wlly,
00793 wurx,wury));
00794 check_nomsg(max_ima1_cy=cpl_image_get_centroid_y_window(img1,wllx,wlly,
00795 wurx,wury));
00796
00797 if(CPL_ERROR_NONE != cpl_image_fit_gaussian(img1,max_ima1_x,max_ima1_y,
00798 SINFO_PSF_SZ,
00799 &norm,&xc,&yc,&sx,&sy,
00800 &fwhm_x,&fwhm_y)) {
00801
00802 sinfo_msg_warning("problem in Gaussian fit");
00803 cpl_error_reset();
00804 }
00805
00806
00807 check_nomsg(bkg1=irplib_strehl_ring_background(img1,max_ima1_x,max_ima1_y,
00808 19,21,
00809 IRPLIB_BG_METHOD_AVER_REJ));
00810
00811
00812 if(CPL_ERROR_NONE != sinfo_strehl_compute(img1,img2,
00813 SINFO_STREHL_M1,SINFO_STREHL_M2,lam,dlam,
00814 pscale1,pscale2,
00815 SINFO_STREHL_BOX_SIZE,
00816 max_ima1_x,max_ima1_y,max_ima2_x,max_ima2_y,
00817 strehl_star_rad1,strehl_bg_rmin1,strehl_bg_rmax1,
00818 NOISE_HSIZE,NOISE_NSAMPLES,
00819 strehl,strehl_err,&star_bkg,&star_peak,&star_flux,
00820 &psf_peak,&psf_flux,&bkg_noise))
00821 {
00822 sinfo_msg_warning("Problem computing strehl, set it to -1");
00823 *strehl=-1;
00824 *strehl_err=0;
00825 irplib_error_reset();
00826 cpl_error_reset();
00827
00828
00829 }
00830
00831 sinfo_free_image(&img1);
00832
00833 return 0;
00834 cleanup:
00835
00836 sinfo_free_image(&img_dup);
00837 return -1;
00838
00839
00840 }
00841
00842
00843
00844
00845 cpl_table* sinfo_get_encircled_energy(cpl_frameset* sof,
00846 cpl_image* img,
00847 double* fwhm_x,
00848 double* fwhm_y,
00849 cpl_table** qclog_tbl)
00850 {
00851
00852 cpl_image* img_dup=NULL;
00853 int max_ima_x=0;
00854 int max_ima_y=0;
00855 int wllx=0;
00856 int wlly=0;
00857 int wurx=0;
00858 int wury=0;
00859 const double d_mirror = 8.;
00860 const double factor = 180/PI_NUMB*3600.;
00861 double max_ima_cx=0;
00862 double max_ima_cy=0;
00863 double norm=0.;
00864 double xc=0.;
00865 double yc=0.;
00866 double sx=0.;
00867 double sy=0.;
00868 double flux=0;
00869 double flux_max=0;
00870 double pix_scale=0;
00871 double lam=0.;
00872 double pscale=0.;
00873 int dr_difr=0;
00874
00875 double r=0.;
00876 double bkg=0.;
00877 int i=0;
00878 int ni=0;
00879 int ir_difr=0;
00880 int dr=0;
00881 int rmin=0;
00882
00883 char band[MAX_NAME_SIZE];
00884 char spat_res[MAX_NAME_SIZE];
00885
00886 cpl_table* enc_energy=NULL;
00887 cpl_frame* frame=NULL;
00888
00889 int ima_szx=0;
00890 int ima_szy=0;
00891
00892
00893
00894 if(NULL != cpl_frameset_find(sof,PRO_COADD_PSF)) {
00895 frame = cpl_frameset_find(sof,PRO_COADD_PSF);
00896 } else if(NULL != cpl_frameset_find(sof,PRO_OBS_PSF)) {
00897 frame = cpl_frameset_find(sof,PRO_OBS_PSF);
00898 } else if(NULL != cpl_frameset_find(sof,PRO_COADD_STD)) {
00899 frame = cpl_frameset_find(sof,PRO_COADD_STD);
00900 } else if(NULL != cpl_frameset_find(sof,PRO_OBS_STD)) {
00901 frame = cpl_frameset_find(sof,PRO_OBS_STD);
00902 } else if(NULL != cpl_frameset_find(sof,PRO_COADD_OBJ)) {
00903 frame = cpl_frameset_find(sof,PRO_COADD_OBJ);
00904 } else if(NULL != cpl_frameset_find(sof,PRO_OBS_OBJ)) {
00905 frame = cpl_frameset_find(sof,PRO_OBS_OBJ);
00906 } else {
00907 sinfo_msg_error("Frame %s or %s or %s or %s or %s or %s not found!",
00908 PRO_COADD_PSF,PRO_OBS_PSF,
00909 PRO_COADD_STD, PRO_OBS_STD,
00910 PRO_COADD_OBJ, PRO_OBS_OBJ);
00911 return NULL;
00912 }
00913
00914 sinfo_get_spatial_res(frame,spat_res);
00915 sinfo_get_band(frame,band);
00916 pix_scale=atof(spat_res);
00917 lam=sinfo_get_wave_cent(band);
00918
00919 pscale=0.5*pix_scale;
00920
00921
00922
00923 dr_difr=factor*1.22*lam*1.e-6/d_mirror/pscale;
00924 ir_difr=floor(dr_difr+0.5);
00925 if (pix_scale==0.025) {
00926 ni=10;
00927 rmin=ir_difr;
00928 dr=rmin;
00929 } else {
00930 ni=15;
00931 sinfo_msg_warning("Reset diffraction limit");
00932 ir_difr=10;
00933 rmin=1;
00934 dr=2;
00935 }
00936
00937 sinfo_msg("Diffraction limit: %d",ir_difr);
00938
00939 check_nomsg(img_dup=cpl_image_duplicate(img));
00940 sinfo_clean_nan(&img_dup);
00941 check_nomsg(cpl_image_get_maxpos(img_dup,&max_ima_x,&max_ima_y));
00942 sinfo_free_image(&img_dup);
00943
00944
00945
00946 check_nomsg(ima_szx=cpl_image_get_size_x(img));
00947 check_nomsg(ima_szy=cpl_image_get_size_y(img));
00948 sinfo_check_borders(&max_ima_x,ima_szx,SINFO_STREHL_WINDOW);
00949 sinfo_check_borders(&max_ima_y,ima_szy,SINFO_STREHL_WINDOW);
00950 sinfo_get_safe_box(&wllx,&wlly,&wurx,&wury,max_ima_x,max_ima_y,SINFO_PSF_SZ,
00951 ima_szx,ima_szy);
00952
00953 check_nomsg(max_ima_cx=cpl_image_get_centroid_x_window(img,wllx,wlly,
00954 wurx,wury));
00955 check_nomsg(max_ima_cy=cpl_image_get_centroid_y_window(img,wllx,wlly,
00956 wurx,wury));
00957
00958 if(CPL_ERROR_NONE != cpl_image_fit_gaussian(img,max_ima_x,max_ima_y,
00959 SINFO_PSF_SZ,
00960 &norm,&xc,&yc,&sx,&sy,
00961 fwhm_x,fwhm_y)) {
00962
00963
00964 sinfo_msg_warning("problem in Gaussian fit");
00965 cpl_error_reset();
00966 }
00967 check_nomsg(enc_energy = cpl_table_new(ni));
00968 check_nomsg(cpl_table_new_column(enc_energy,"r_pix", CPL_TYPE_INT));
00969 check_nomsg(cpl_table_new_column(enc_energy,"r_mas", CPL_TYPE_DOUBLE));
00970 check_nomsg(cpl_table_new_column(enc_energy,"r_dif", CPL_TYPE_DOUBLE));
00971 check_nomsg(cpl_table_new_column(enc_energy,"abs_energy" , CPL_TYPE_DOUBLE));
00972 check_nomsg(cpl_table_new_column(enc_energy,"rel_energy" , CPL_TYPE_DOUBLE));
00973
00974 check_nomsg(bkg=irplib_strehl_ring_background(img,max_ima_x,max_ima_y,
00975 25,30,IRPLIB_BG_METHOD_AVER_REJ)) ;
00976 r=rmin+(ni-1)*dr;
00977 check_nomsg(flux_max=irplib_strehl_disk_flux(img,max_ima_x,max_ima_y,r,bkg));
00978 r=rmin;
00979
00980 for(i=0; i<ni; i++)
00981 {
00982 check_nomsg(flux=irplib_strehl_disk_flux(img,max_ima_x,max_ima_y,r,bkg));
00983 check_nomsg(cpl_table_set_int(enc_energy,"r_pix",i,r));
00984 check_nomsg(cpl_table_set_double(enc_energy,"r_mas",i,r*pscale));
00985 check_nomsg(cpl_table_set_double(enc_energy,"r_dif",i,r/ir_difr));
00986 check_nomsg(cpl_table_set_double(enc_energy,"abs_energy",i,flux));
00987 check_nomsg(cpl_table_set_double(enc_energy,"rel_energy",i,flux/flux_max));
00988 r+=dr;
00989
00990 }
00991
00992 sinfo_msg("max ima=%d %d\n",max_ima_x,max_ima_y);
00993 sinfo_msg("centroid ima=%f %f\n",max_ima_cx,max_ima_cy);
00994 sinfo_msg("gauss info=%f %f %f %f %f %f %f\n",
00995 norm,xc,yc,sx,sy,*fwhm_x,*fwhm_y);
00996
00997 check_nomsg(flux=irplib_strehl_disk_flux(img,max_ima_x,max_ima_y,
00998 ir_difr,bkg));
00999 ck0_nomsg(sinfo_qclog_add_double(*qclog_tbl,"QC ENC CORE",
01000 flux/flux_max,
01001 "Encircled energy within PSF core","%f"));
01002
01003 return enc_energy;
01004
01005 cleanup:
01006 sinfo_free_image(&img_dup);
01007
01008 return NULL;
01009 }
01010
01011
01012 static cpl_table* sinfo_get_strehl_from_cube(cpl_imagelist* cube,
01013 char* name,
01014 cpl_frame* frame)
01015 {
01016 cpl_table* strehl_tbl=NULL;
01017
01018 double dispersion=0.;
01019 double centralWave=0.;
01020 double wrange=0;
01021 double wstart=0;
01022 double wstep=0;
01023 double wend=0;
01024 double ws=0;
01025 double we=0;
01026 double pix_scale=0;
01027 double lam=0;
01028 double dlam=0;
01029 double pscale = 0;
01030
01031 double strehl_star_radius=0;
01032 double strehl_bg_r1=0;
01033 double strehl_bg_r2=0;
01034 double strehl=0;
01035 double strehl_err=0;
01036 char spat_res[MAX_NAME_SIZE];
01037 cpl_propertylist* plist=NULL;
01038
01039 int naxis3=0;
01040 int nsample=0;
01041 int i=0;
01042
01043
01044 sinfo_get_spatial_res(frame,spat_res);
01045 pix_scale=atof(spat_res);
01046 sinfo_msg("pix scale=%f",pix_scale);
01047
01048 pscale=0.5*pix_scale;
01049
01050 strehl_star_radius=25*pscale;
01051 strehl_bg_r1=25*pscale;
01052 strehl_bg_r2=30*pscale;
01053
01054 plist=cpl_propertylist_load(name,0);
01055 dispersion=sinfo_pfits_get_cdelt3(plist);
01056 centralWave=sinfo_pfits_get_crval3(plist);
01057 naxis3=sinfo_pfits_get_naxis3(plist);
01058 sinfo_free_propertylist(&plist);
01059 wrange=dispersion*naxis3;
01060
01061 wstart = centralWave - (float) (cpl_imagelist_get_size(cube) / 2)*
01062 dispersion+dispersion;
01063 wend =wstart + dispersion * cpl_imagelist_get_size(cube);
01064 wstep=0.025;
01065
01066
01067
01068
01069
01070 nsample=(int)((wend-wstart-wstep)/wstep);
01071 check_nomsg(strehl_tbl = cpl_table_new(nsample));
01072 check_nomsg(cpl_table_new_column(strehl_tbl,"wavelength",CPL_TYPE_DOUBLE));
01073 check_nomsg(cpl_table_new_column(strehl_tbl,"strehl",CPL_TYPE_DOUBLE));
01074 check_nomsg(cpl_table_new_column(strehl_tbl,"strehl_error",CPL_TYPE_DOUBLE));
01075
01076
01077 for(i=1;i<nsample;i++) {
01078
01079 ws=wstart+wstep*i;
01080 we=ws+wstep;
01081
01082 lam = (double)0.5*(ws+we);
01083 dlam=wstep;
01084
01085 check(sinfo_get_strehl_from_slice(cube,
01086 dispersion,
01087 centralWave,
01088 ws,
01089 we,
01090 pscale,
01091 strehl_star_radius,
01092 strehl_bg_r1,
01093 strehl_bg_r2,
01094 &strehl,
01095 &strehl_err),"Error computing strehl");
01096
01097 if((isnan(lam) ==0) &&
01098 (isnan(lam) ==0) &&
01099 (isnan(lam) ==0)) {
01100 check_nomsg(cpl_table_set_double(strehl_tbl,"wavelength",i,lam));
01101 check_nomsg(cpl_table_set_double(strehl_tbl,"strehl",i,strehl));
01102 check_nomsg(cpl_table_set_double(strehl_tbl,"strehl_error",i,
01103 strehl_err));
01104
01105 }
01106 }
01107
01108 return strehl_tbl;
01109
01110 cleanup:
01111 return NULL;
01112
01113
01114 }
01115
01116 static double sinfo_get_strehl_from_ima(cpl_image* ima,
01117 char* name,
01118 cpl_frame* frame)
01119 {
01120
01121 double dispersion=0.;
01122 double centralWave=0.;
01123 double wrange=0;
01124 double wstart=0;
01125 double wend=0;
01126 double pix_scale=0;
01127 double pscale = 0;
01128
01129 double strehl_star_radius=0;
01130 double strehl_bg_r1=0;
01131 double strehl_bg_r2=0;
01132 double strehl=0;
01133 double strehl_err=0;
01134 char spat_res[MAX_NAME_SIZE];
01135
01136 int naxis3=0;
01137 cpl_propertylist* plist=NULL;
01138
01139
01140 sinfo_get_spatial_res(frame,spat_res);
01141 pix_scale=atof(spat_res);
01142
01143 pscale=0.5*pix_scale;
01144
01145
01146
01147 strehl_star_radius=25*pscale;
01148 strehl_bg_r1=25*pscale;
01149 strehl_bg_r2=30*pscale;
01150
01151 plist=cpl_propertylist_load(name,0);
01152 dispersion=sinfo_pfits_get_cdelt3(plist);
01153 centralWave=sinfo_pfits_get_crval3(plist);
01154 naxis3=sinfo_pfits_get_naxis3(plist);
01155 sinfo_free_propertylist(&plist);
01156 wrange=dispersion*naxis3;
01157
01158 wstart = centralWave - (wrange / 2);
01159 wend = centralWave + (wrange / 2);
01160
01161 check(sinfo_get_strehl_from_image(ima,
01162 wstart,
01163 wend,
01164 pscale,
01165 strehl_star_radius,
01166 strehl_bg_r1,
01167 strehl_bg_r2,
01168 &strehl,
01169 &strehl_err),"Computing Strehl");
01170
01171
01172
01173
01174
01175 cleanup:
01176 return strehl;
01177
01178
01179 }
01180
01181 static int
01182 sinfo_get_frm12(cpl_frameset* sof,cpl_frame** frm1,cpl_frame** frm2){
01183
01184 cpl_frameset* obs=NULL;
01185 int nobs=0;
01186 float eps=0.0001;
01187 float* pix_scale=NULL;
01188 int i=0;
01189 cpl_frame* frame=NULL;
01190
01191 obs = cpl_frameset_new();
01192 sinfo_contains_frames_kind(sof,obs,PRO_OBS_PSF);
01193 nobs=cpl_frameset_get_size(obs);
01194 if (nobs < 1) {
01195 sinfo_contains_frames_kind(sof,obs,PRO_OBS_STD);
01196 nobs=cpl_frameset_get_size(obs);
01197 }
01198
01199 nobs=cpl_frameset_get_size(obs);
01200 if (nobs < 1) {
01201 return -1;
01202 } else {
01203 pix_scale=cpl_calloc(nobs,sizeof(float));
01204 for(i=0;i<nobs;i++) {
01205 frame=cpl_frameset_get_frame(obs,i);
01206 pix_scale[i]=sinfo_pfits_get_pixelscale(
01207 (char*)cpl_frame_get_filename(frame));
01208 if(fabs(pix_scale[i]-0.025)< eps) {
01209 *frm1=cpl_frame_duplicate(frame);
01210 } else if (fabs(pix_scale[i]-0.1) <eps) {
01211 *frm2=cpl_frame_duplicate(frame);
01212 } else {
01213 sinfo_msg_error("No proper frame found for strehl computation");
01214 return -1;
01215 }
01216 }
01217 }
01218 cpl_free(pix_scale);
01219 cpl_frameset_delete(obs);
01220
01221 return 0;
01222
01223 }
01224
01225 static cpl_table*
01226 sinfo_get_strehl_from_2cubes(cpl_imagelist* cube1,
01227 cpl_imagelist* cube2,
01228 const char* fname1,
01229 cpl_frame* frm1,
01230 cpl_frame* frm2)
01231 {
01232
01233 cpl_table* strehl_tbl=NULL;
01234
01235
01236 double dispersion=0.;
01237 double centralWave=0.;
01238 double wrange=0;
01239 double wstart=0;
01240 double wstep=0;
01241 double wend=0;
01242 double ws=0;
01243 double we=0;
01244 double pix_scale1=0;
01245 double pix_scale2=0;
01246 double lam=0;
01247 double dlam=0;
01248 double pscale1 = 0;
01249 double pscale2 = 0;
01250
01251 double strehl_star_rad1=0;
01252 double strehl_star_rad2=0;
01253 double strehl_bg_rmin1=0;
01254 double strehl_bg_rmin2=0;
01255 double strehl_bg_rmax1=0;
01256 double strehl_bg_rmax2=0;
01257 double strehl=0;
01258 double strehl_err=0;
01259 char res1[MAX_NAME_SIZE];
01260 char res2[MAX_NAME_SIZE];
01261
01262 int naxis3=0;
01263 int nsample=0;
01264 int i=0;
01265 cpl_propertylist* plist=NULL;
01266
01267 sinfo_get_spatial_res(frm1,res1);
01268 sinfo_get_spatial_res(frm2,res2);
01269 pix_scale1=atof(res1);
01270 pix_scale2=atof(res2);
01271
01272 pscale1=0.5*pix_scale1;
01273 pscale2=0.5*pix_scale2;
01274
01275
01276
01277 strehl_star_rad1=25*pscale1;
01278 strehl_bg_rmin1=25*pscale1;
01279 strehl_bg_rmax1=30*pscale1;
01280
01281 strehl_star_rad2=25*pscale2;
01282 strehl_bg_rmin2=25*pscale2;
01283 strehl_bg_rmax2=30*pscale2;
01284 plist=cpl_propertylist_load(fname1,0);
01285 dispersion=sinfo_pfits_get_cdelt3(plist);
01286 centralWave=sinfo_pfits_get_crval3(plist);
01287 naxis3=sinfo_pfits_get_naxis3(plist);
01288 sinfo_free_propertylist(&plist);
01289 wrange=dispersion*naxis3;
01290
01291 wstart = centralWave - (float) (cpl_imagelist_get_size(cube1) / 2)*
01292 dispersion+dispersion;
01293 wend = wstart + dispersion * cpl_imagelist_get_size(cube1);
01294 wstep = 0.025;
01295
01296
01297
01298
01299
01300 nsample=(int)((wend-wstart-wstep)/wstep);
01301
01302 check_nomsg(strehl_tbl = cpl_table_new(nsample));
01303 check_nomsg(cpl_table_new_column(strehl_tbl,"wavelength",CPL_TYPE_DOUBLE));
01304 check_nomsg(cpl_table_new_column(strehl_tbl,"strehl",CPL_TYPE_DOUBLE));
01305 check_nomsg(cpl_table_new_column(strehl_tbl,"strehl_error",CPL_TYPE_DOUBLE));
01306
01307
01308 for(i=1;i<nsample;i++) {
01309
01310 ws=wstart+wstep*i;
01311 we=ws+wstep;
01312
01313 lam = (double)0.5*(ws+we);
01314 dlam=wstep;
01315
01316 sinfo_get_strehl_from_2slices(cube1,
01317 cube2,
01318 dispersion,
01319 centralWave,
01320 ws,
01321 we,
01322 pscale1,
01323 pscale2,
01324 strehl_star_rad1,
01325 strehl_bg_rmin1,
01326 strehl_bg_rmax1,
01327 &strehl,
01328 &strehl_err);
01329
01330
01331 if((isnan(lam) ==0) &&
01332 (isnan(lam) ==0) &&
01333 (isnan(lam) ==0)) {
01334 check_nomsg(cpl_table_set_double(strehl_tbl,"wavelength",i,lam));
01335 check_nomsg(cpl_table_set_double(strehl_tbl,"strehl",i,strehl));
01336 check_nomsg(cpl_table_set_double(strehl_tbl,"strehl_error",
01337 i,strehl_err));
01338
01339 }
01340
01341 }
01342
01343 return strehl_tbl;
01344 cleanup:
01345 return NULL;
01346 }
01347
01348
01349
01350
01383
01384 int sinfo_strehl_compute(
01385 const cpl_image * im1,
01386 const cpl_image * im2,
01387 double m1,
01388 double m2,
01389 double lam,
01390 double dlam,
01391 double pscale1,
01392 double pscale2,
01393 int size,
01394 int xpos1,
01395 int ypos1,
01396 int xpos2,
01397 int ypos2,
01398 double r1,
01399 double r2,
01400 double r3,
01401 int noise_box_sz,
01402 int noise_nsamples,
01403 double * strehl,
01404 double * strehl_err,
01405 double * star_bg,
01406 double * star_peak,
01407 double * star_flux,
01408 double * psf_peak,
01409 double * psf_flux,
01410 double * bg_noise)
01411 {
01412 cpl_image * psf ;
01413 const int window_size = 5;
01414 int star_radius ;
01415 int ring[4] ;
01416 double peak2 =0;
01417
01418
01419 irplib_assure_code(im1 != NULL, CPL_ERROR_NULL_INPUT);
01420 irplib_assure_code(im2 != NULL, CPL_ERROR_NULL_INPUT);
01421 irplib_assure_code(strehl != NULL, CPL_ERROR_NULL_INPUT);
01422 irplib_assure_code(strehl_err != NULL, CPL_ERROR_NULL_INPUT);
01423 irplib_assure_code(star_bg != NULL, CPL_ERROR_NULL_INPUT);
01424 irplib_assure_code(star_peak != NULL, CPL_ERROR_NULL_INPUT);
01425 irplib_assure_code(star_flux != NULL, CPL_ERROR_NULL_INPUT);
01426 irplib_assure_code(psf_peak != NULL, CPL_ERROR_NULL_INPUT);
01427 irplib_assure_code(psf_flux != NULL, CPL_ERROR_NULL_INPUT);
01428
01429 irplib_assure_code(pscale1 > 0.0, CPL_ERROR_ILLEGAL_INPUT);
01430 irplib_assure_code(pscale2 > 0.0, CPL_ERROR_ILLEGAL_INPUT);
01431
01432 irplib_assure_code(xpos1-window_size > 0, CPL_ERROR_ACCESS_OUT_OF_RANGE);
01433 irplib_assure_code(ypos1-window_size > 0, CPL_ERROR_ACCESS_OUT_OF_RANGE);
01434 irplib_assure_code(xpos2-window_size > 0, CPL_ERROR_ACCESS_OUT_OF_RANGE);
01435 irplib_assure_code(ypos2-window_size > 0, CPL_ERROR_ACCESS_OUT_OF_RANGE);
01436
01437 irplib_assure_code(xpos1+window_size <= cpl_image_get_size_x(im1),
01438 CPL_ERROR_ACCESS_OUT_OF_RANGE);
01439 irplib_assure_code(ypos1+window_size <= cpl_image_get_size_y(im1),
01440 CPL_ERROR_ACCESS_OUT_OF_RANGE);
01441
01442 irplib_assure_code(xpos2+window_size <= cpl_image_get_size_x(im2),
01443 CPL_ERROR_ACCESS_OUT_OF_RANGE);
01444 irplib_assure_code(ypos2+window_size <= cpl_image_get_size_y(im2),
01445 CPL_ERROR_ACCESS_OUT_OF_RANGE);
01446
01447 irplib_assure_code(r1 > 0.0, CPL_ERROR_ILLEGAL_INPUT);
01448 irplib_assure_code(r2 > 0.0, CPL_ERROR_ILLEGAL_INPUT);
01449 irplib_assure_code(r3 > r2, CPL_ERROR_ILLEGAL_INPUT);
01450
01451
01452
01453
01454
01455
01456 psf = irplib_strehl_generate_psf(m1, m2, lam, dlam, pscale2, size);
01457 irplib_assure_code(psf != NULL, CPL_ERROR_ILLEGAL_OUTPUT);
01458
01459
01460 *psf_peak = cpl_image_get_max(psf) ;
01461
01462
01463
01464 cpl_image_delete(psf) ;
01465 assert( *psf_peak > 0.0);
01466 *psf_flux = 1.0;
01467
01468 ring[0] = xpos2 ;
01469 ring[1] = ypos2 ;
01470 ring[2] = (int)(r2/pscale2);
01471 ring[3] = (int)(r3/pscale2);
01472
01473
01474 *star_bg = irplib_strehl_ring_background(im2, xpos2, ypos2,
01475 ring[2], ring[3],
01476 IRPLIB_BG_METHOD_AVER_REJ);
01477
01478
01479 star_radius = (int)(r1/pscale2);
01480
01481
01482 *star_flux = irplib_strehl_disk_flux(im2, xpos2, ypos2, star_radius,
01483 *star_bg);
01484
01485 irplib_assure_code(*star_flux > 0.0, CPL_ERROR_ILLEGAL_OUTPUT);
01486
01487
01488
01489
01490 *star_peak = cpl_image_get_max_window(im1,
01491 xpos1-window_size,
01492 ypos1-window_size,
01493 xpos1+window_size,
01494 ypos1+window_size) - *star_bg ;
01495 peak2 = cpl_image_get_max_window(im2,
01496 xpos2-window_size,
01497 ypos2-window_size,
01498 xpos2+window_size,
01499 ypos2+window_size) - *star_bg ;
01500 if ( fabs(peak2-*star_peak)/peak2 > 0.25) {
01501 sinfo_msg_warning("rel diff: %g",fabs(peak2-*star_peak)/peak2);
01502 }
01503 irplib_assure_code(*star_peak > 0.0, CPL_ERROR_ILLEGAL_OUTPUT);
01504
01505
01506 *strehl = (*star_peak / *star_flux) / (*psf_peak / *psf_flux) ;
01507
01508 if (*strehl > 1)
01509 cpl_msg_warning(cpl_func, "Extreme Strehl-ratio=%g, star_peak=%g, "
01510 "star_flux=%g, psf_peak=%g, psf_flux=%g", *strehl,
01511 *star_peak, *star_flux, *psf_peak, *psf_flux);
01512
01513
01514 if (cpl_flux_get_noise_ring(im2, ring, noise_box_sz, noise_nsamples,
01515 bg_noise, NULL) == CPL_ERROR_NONE) {
01516 *strehl_err = SINFO_STREHL_ERROR_COEFFICIENT * (*bg_noise) * pscale2 *
01517 star_radius * star_radius / *star_flux;
01518 irplib_assure_code(*strehl_err >= 0.0, CPL_ERROR_ILLEGAL_OUTPUT);
01519 }
01520
01521 return CPL_ERROR_NONE;
01522 }
01523
01524 static void
01525 sinfo_check_borders(int* val,const int max,const int thresh)
01526 {
01527
01528 *val = ((*val-thresh) > 0) ? *val : thresh;
01529 *val = ((*val+thresh) < max) ? *val : max-thresh-1;
01530 return;
01531 }
01532
01533 static void
01534 sinfo_get_safe_box(int* llx,
01535 int* lly,
01536 int* urx,
01537 int* ury,
01538 const int xpos,
01539 const int ypos,
01540 const int box,
01541 const int szx,
01542 const int szy)
01543
01544 {
01545 *llx= ((xpos-box)>0) ? (xpos-box) : 1;
01546 *lly= ((ypos-box)>0) ? (ypos-box) : 1;
01547 *urx= ((xpos+box)<szx) ? (xpos+box) : szx-1 ;
01548 *ury= ((ypos+box)<szy) ? (ypos+box) : szy-1 ;
01549 return;
01550 }