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 #include <fors_extract.h>
00033
00034 #include <fors_tools.h>
00035 #include <fors_dfs.h>
00036 #include <fors_pfits.h>
00037 #include <fors_utils.h>
00038
00039 #include <cpl.h>
00040
00041 #include <string.h>
00042 #include <stdbool.h>
00043
00050 struct _extract_method
00051 {
00052 enum {SEX, TEST} method;
00053 const char *sex_exe;
00054 const char *sex_config;
00055 const char *sex_mag;
00056 const char *sex_magerr;
00057 int sex_radius;
00058 };
00059
00060 static fors_star_list *
00061 extract_sex(const fors_image *image,
00062 const fors_setting *setting,
00063 const char *sex_exe,
00064 const char *sex_config,
00065 const char *sex_mag,
00066 const char *sex_magerr,
00067 int radius,
00068 fors_extract_sky_stats *sky_stats,
00069 cpl_image **background,
00070 cpl_table **extracted_sources);
00071
00072 static fors_star_list *
00073 extract_test(fors_extract_sky_stats *sky_stats,
00074 cpl_image **background,
00075 cpl_table **extracted_sources);
00076
00082 void
00083 fors_extract_define_parameters(cpl_parameterlist *parameters,
00084 const char *context)
00085 {
00086 cpl_parameter *p;
00087 const char *full_name = NULL;
00088 const char *name;
00089
00090 name = "extract_method";
00091 full_name = cpl_sprintf("%s.%s", context, name);
00092 p = cpl_parameter_new_enum(full_name,
00093 CPL_TYPE_STRING,
00094 "Source extraction method",
00095 context,
00096 "sex", 2,
00097 "sex", "test");
00098 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00099 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00100 cpl_parameterlist_append(parameters, p);
00101 cpl_free((void *)full_name);
00102
00103 name = "sex_exe";
00104 full_name = cpl_sprintf("%s.%s", context, name);
00105 p = cpl_parameter_new_value(full_name,
00106 CPL_TYPE_STRING,
00107 "SExtractor executable",
00108 context,
00109 FORS_SEXTRACTOR_PATH "/sex");
00110 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00111 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00112 cpl_parameterlist_append(parameters, p);
00113 cpl_free((void *)full_name);
00114
00115 name = "sex_config";
00116 full_name = cpl_sprintf("%s.%s", context, name);
00117 p = cpl_parameter_new_value(full_name,
00118 CPL_TYPE_STRING,
00119 "SExtractor configuration file",
00120 context,
00121 FORS_SEXTRACTOR_CONFIG "/fors.sex");
00122 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00123 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00124 cpl_parameterlist_append(parameters, p);
00125 cpl_free((void *)full_name);
00126
00127
00128 name = "sex_mag";
00129 full_name = cpl_sprintf("%s.%s", context, name);
00130 p = cpl_parameter_new_value(full_name,
00131 CPL_TYPE_STRING,
00132 "SExtractor magnitude",
00133 context,
00134 "MAG_APER");
00135 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00136 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00137 cpl_parameterlist_append(parameters, p);
00138 cpl_free((void *)full_name);
00139
00140 name = "sex_magerr";
00141 full_name = cpl_sprintf("%s.%s", context, name);
00142 p = cpl_parameter_new_value(full_name,
00143 CPL_TYPE_STRING,
00144 "SExtractor magnitude error",
00145 context,
00146 "MAGERR_APER");
00147 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00148 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00149 cpl_parameterlist_append(parameters, p);
00150 cpl_free((void *)full_name);
00151
00152 name = "sex_radius";
00153 full_name = cpl_sprintf("%s.%s", context, name);
00154 p = cpl_parameter_new_value(full_name,
00155 CPL_TYPE_INT,
00156 "Background error map median filter "
00157 "radius (unbinned pixels)",
00158 context,
00159 64);
00160 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00161 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00162 cpl_parameterlist_append(parameters, p);
00163 cpl_free((void *)full_name); full_name = NULL;
00164
00165 return;
00166 }
00167
00168 #undef cleanup
00169 #define cleanup \
00170 do { \
00171 cpl_free((void *)name); \
00172 } while (0)
00173
00182 extract_method *
00183 fors_extract_method_new(const cpl_parameterlist *parameters, const char *context)
00184 {
00185 extract_method *em = cpl_malloc(sizeof(*em));
00186 const char *name = NULL;
00187 const char *method;
00188
00189 cpl_msg_info(cpl_func, "Extraction method:");
00190
00191 cpl_msg_indent_more();
00192 name = cpl_sprintf("%s.%s", context, "extract_method");
00193 method = dfs_get_parameter_string_const(parameters,
00194 name);
00195 cpl_msg_indent_less();
00196
00197 assure( !cpl_error_get_code(), return NULL, NULL );
00198 assure( method != NULL, return NULL, NULL );
00199
00200 if (strcmp(method, "sex") == 0) {
00201 em->method = SEX;
00202
00203 cpl_msg_indent_more();
00204 name = cpl_sprintf("%s.%s", context, "sex_exe");
00205 em->sex_exe = dfs_get_parameter_string_const(parameters,
00206 name);
00207 cpl_free((void *)name); name = NULL;
00208 cpl_msg_indent_less();
00209
00210
00211 cpl_msg_indent_more();
00212 name = cpl_sprintf("%s.%s", context, "sex_config");
00213 em->sex_config = dfs_get_parameter_string_const(parameters,
00214 name);
00215 cpl_free((void *)name); name = NULL;
00216 cpl_msg_indent_less();
00217
00218
00219
00220 cpl_msg_indent_more();
00221 name = cpl_sprintf("%s.%s", context, "sex_mag");
00222 em->sex_mag = dfs_get_parameter_string_const(parameters,
00223 name);
00224 cpl_free((void *)name); name = NULL;
00225 cpl_msg_indent_less();
00226
00227
00228 cpl_msg_indent_more();
00229 name = cpl_sprintf("%s.%s", context, "sex_magerr");
00230 em->sex_magerr = dfs_get_parameter_string_const(parameters,
00231 name);
00232 cpl_free((void *)name); name = NULL;
00233 cpl_msg_indent_less();
00234
00235
00236 cpl_msg_indent_more();
00237 name = cpl_sprintf("%s.%s", context, "sex_radius");
00238 em->sex_radius = dfs_get_parameter_int_const(parameters,
00239 name);
00240 cpl_free((void *)name); name = NULL;
00241 cpl_msg_indent_less();
00242 }
00243 else if (strcmp(method, "test") == 0) {
00244 em->method = TEST;
00245 }
00246 else {
00247 assure( false, return NULL, "Unknown extraction method '%s'", method);
00248 }
00249
00250 cleanup;
00251 return em;
00252 }
00253
00257 void
00258 fors_extract_method_delete(extract_method **em)
00259 {
00260 if (em && *em) {
00261 cpl_free(*em); *em = NULL;
00262 }
00263 return;
00264 }
00265
00266 #undef cleanup
00267 #define cleanup
00268
00278 fors_star_list *
00279 fors_extract(const fors_image *image,
00280 const fors_setting *setting,
00281 const extract_method *em,
00282 fors_extract_sky_stats *sky_stats,
00283 cpl_image **background,
00284 cpl_table **extracted_sources)
00285 {
00286 assure( em != NULL, return NULL, NULL );
00287
00288 cpl_msg_info(cpl_func, "Extracting sources");
00289
00290 switch (em->method ) {
00291 case SEX: return extract_sex(image, setting,
00292 em->sex_exe,
00293 em->sex_config,
00294 em->sex_mag,
00295 em->sex_magerr,
00296 em->sex_radius,
00297 sky_stats,
00298 background,
00299 extracted_sources); break;
00300 case TEST: return extract_test(sky_stats, background, extracted_sources); break;
00301 default:
00302 assure( false, return NULL, "Unknown method %d", em->method );
00303 break;
00304 }
00305 }
00306
00307 #undef cleanup
00308 #define cleanup \
00309 do { \
00310 cpl_table_delete(out); \
00311 cpl_free((void *)command); \
00312 } while (0)
00313
00337 static fors_star_list *
00338 extract_sex(const fors_image *image,
00339 const fors_setting *setting,
00340 const char *sex_exe,
00341 const char *sex_config,
00342 const char *sex_mag,
00343 const char *sex_magerr,
00344 int radius,
00345 fors_extract_sky_stats *sky_stats,
00346 cpl_image **background,
00347 cpl_table **extracted_sources)
00348 {
00349 const char *const filename_data = "sextract_data.fits";
00350 const char *const filename_sigma = "sextract_bkg_sigma.fits";
00351 const char *const filename_cat = "sextract_cat.fits";
00352 const char *const filename_bkg = "sextract_bkg.fits";
00353 cpl_table *out = NULL;
00354 const char *command = NULL;
00355 fors_star_list *stars = NULL;
00356
00357 assure( setting != NULL, return NULL, NULL );
00358
00359 assure( image != NULL, return NULL, NULL );
00360
00361 assure( sky_stats != NULL, return NULL, NULL );
00362 assure( background != NULL, return NULL, NULL );
00363
00364
00365
00366 fors_image_save_sex(image, NULL,
00367 filename_data,
00368 filename_sigma,
00369 radius);
00370 assure( !cpl_error_get_code(), return NULL,
00371 "Could not save image to %s and %s",
00372 filename_data, filename_sigma);
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383 command = cpl_sprintf("%s %s,%s "
00384 "-c %s "
00385
00386
00387
00388 "-GAIN %f "
00389
00390
00391
00392 "-CHECKIMAGE_TYPE BACKGROUND "
00393 "-CHECKIMAGE_NAME %s "
00394 "-WEIGHT_TYPE MAP_RMS "
00395 "-WEIGHT_IMAGE %s,%s "
00396 "-CATALOG_TYPE FITS_1.0 "
00397 "-CATALOG_NAME %s",
00398 sex_exe,
00399 filename_data, filename_data,
00400 sex_config,
00401 1.0/setting->average_gain,
00402 filename_bkg,
00403 filename_sigma, filename_sigma,
00404 filename_cat);
00405
00406 cpl_msg_info(cpl_func, "Running '%s'", command);
00407
00408 assure( system(command) == 0, return stars, "'%s' failed", command);
00409
00410
00411
00412
00413
00414 {
00415 int plane = 0;
00416 int extension = 0;
00417 cpl_image *work_back;
00418
00419 *background = cpl_image_load(filename_bkg,
00420 CPL_TYPE_FLOAT, plane, extension);
00421 assure( !cpl_error_get_code(), return NULL,
00422 "Could not load SExtractor background image %s",
00423 filename_bkg );
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434 int nx = cpl_image_get_size_x(*background);
00435 int ny = cpl_image_get_size_y(*background);
00436 int xlo = nx/2 - 50;
00437 int xhi = nx/2 + 50;
00438 int ylo = ny/2 - 50;
00439 int yhi = ny/2 + 50;
00440
00441 if (xlo < 0) xlo = 0;
00442 if (xhi >= nx) xhi = nx - 1;
00443 if (ylo < 0) ylo = 0;
00444 if (yhi >= ny) yhi = ny - 1;
00445
00446 work_back = cpl_image_duplicate(*background);
00447
00448 sky_stats->mean = cpl_image_get_mean_window(work_back,
00449 xlo, ylo, xhi, yhi);
00450 sky_stats->median = cpl_image_get_median_window(work_back,
00451 xlo, ylo, xhi, yhi);
00452 cpl_image_subtract_scalar(work_back, sky_stats->median);
00453 cpl_image_abs(work_back);
00454 sky_stats->rms = cpl_image_get_median_window(work_back,
00455 xlo, ylo, xhi, yhi)
00456 * STDEV_PR_MAD;
00457
00458 cpl_image_delete(work_back);
00459
00460 assure( !cpl_error_get_code(), return NULL,
00461 "Could not calculate sky statistics" );
00462
00463 }
00464
00465 cpl_msg_info(cpl_func, "Background = %f +- %f ADU",
00466 sky_stats->median, sky_stats->rms);
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480 float level;
00481 cpl_image *bmaxsigma;
00482 {
00483 int xradius = 64;
00484 int yradius = 64;
00485 int plane = 0;
00486 int extension = 0;
00487 cpl_image *bsigma;
00488 fors_image *fbsigma;
00489 float maxima;
00490 float minima;
00491
00492 bsigma = cpl_image_load(filename_sigma,
00493 CPL_TYPE_FLOAT, plane, extension);
00494
00495 assure( !cpl_error_get_code(), return NULL,
00496 "Could not load SExtractor background error image %s",
00497 filename_bkg );
00498
00499
00500
00501
00502
00503
00504 fbsigma = fors_image_new(cpl_image_duplicate(bsigma), bsigma);
00505
00506 {
00507 bool use_variance = true;
00508 bmaxsigma = fors_image_filter_max_create(fbsigma, xradius,
00509 yradius, use_variance);
00510
00511
00512
00513
00514 }
00515
00516 fors_image_delete(&fbsigma);
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527 level = cpl_image_get_median(bmaxsigma) * 5;
00528
00529 cpl_msg_debug(cpl_func, "Threshold level = %f",
00530 level);
00531
00532 }
00533
00534 out = cpl_table_load(filename_cat, 1, 1);
00535
00536 assure( !cpl_error_get_code(), return NULL,
00537 "Could not load SExtractor output table %s",
00538 filename_cat);
00539
00540
00541 assure( cpl_table_has_column(out, "FLAGS"), return NULL,
00542 "%s: Missing column: %s", filename_cat, "FLAGS");
00543
00544 assure( cpl_table_has_column(out, "CLASS_STAR"), return NULL,
00545 "%s: Missing column: %s", filename_cat, "CLASS_STAR");
00546
00547 assure( cpl_table_has_column(out, "BACKGROUND"), return NULL,
00548 "%s: Missing column: %s", filename_cat, "BACKGROUND");
00549
00550 assure( cpl_table_has_column(out, "X_IMAGE"), return NULL,
00551 "%s: Missing column: %s", filename_cat, "X_IMAGE");
00552
00553 assure( cpl_table_has_column(out, "Y_IMAGE"), return NULL,
00554 "%s: Missing column: %s", filename_cat, "Y_IMAGE");
00555
00556 assure( cpl_table_has_column(out, "FWHM_IMAGE"), return NULL,
00557 "%s: Missing column: %s", filename_cat, "FWHM_IMAGE");
00558
00559 assure( cpl_table_has_column(out, "A_IMAGE"), return NULL,
00560 "%s: Missing column: %s", filename_cat, "A_IMAGE");
00561
00562 assure( cpl_table_has_column(out, "B_IMAGE"), return NULL,
00563 "%s: Missing column: %s", filename_cat, "B_IMAGE");
00564
00565 assure( cpl_table_has_column(out, "THETA_IMAGE"), return NULL,
00566 "%s: Missing column: %s", filename_cat, "THETA_IMAGE");
00567
00568 assure( cpl_table_has_column(out, sex_mag), return NULL,
00569 "%s: Missing column: %s", filename_cat, sex_mag);
00570
00571 assure( cpl_table_has_column(out, sex_magerr), return NULL,
00572 "%s: Missing column: %s", filename_cat, sex_magerr);
00573
00574
00575
00576
00577 stars = fors_star_list_new();
00578
00579 {
00580 int i;
00581 int bkg_rejected = 0;
00582 int rejected = 0;
00583 float *bdata = cpl_image_get_data(bmaxsigma);
00584 int nx = cpl_image_get_size_x(bmaxsigma);
00585 int ny = cpl_image_get_size_y(bmaxsigma);
00586
00587 for (i = 0; i < cpl_table_get_nrow(out); i++) {
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599 int xpos = cpl_table_get_float(out, "X_IMAGE", i, NULL);
00600 int ypos = cpl_table_get_float(out, "Y_IMAGE", i, NULL);
00601
00602 int flag = cpl_table_get_int(out, "FLAGS", i, NULL);
00603 float fwhm = cpl_table_get_float(out, "FWHM_IMAGE", i, NULL);
00604 float magnitude = cpl_table_get_float(out, sex_mag, i, NULL);
00605
00606 cpl_msg_debug(cpl_func,
00607 "Source at (%d, %d): flag = %d; fwhm = %f pix; "
00608 "background error = %f; level = %f",
00609 xpos, ypos, flag, fwhm,
00610 (xpos > 0 && xpos <= nx && ypos > 0 && ypos <= ny) ?
00611 bdata[(xpos-1) + (ypos-1)*nx] :
00612 -1, level);
00613
00614 if (flag == 0 && fwhm >= 0 && magnitude < 98.0) {
00615
00616 if (xpos > 0 && xpos <= nx && ypos > 0 && ypos <= ny &&
00617 bdata[(xpos-1) + (ypos-1)*nx] < level) {
00618
00619
00620
00621
00622
00623 fors_star *s = fors_star_new(
00624 cpl_table_get_float(out, "X_IMAGE", i, NULL),
00625 cpl_table_get_float(out, "Y_IMAGE", i, NULL),
00626 cpl_table_get_float(out, "FWHM_IMAGE", i, NULL),
00627 cpl_table_get_float(out, "A_IMAGE", i, NULL),
00628 cpl_table_get_float(out, "B_IMAGE", i, NULL),
00629 cpl_table_get_float(out, "THETA_IMAGE",
00630 i, NULL) * M_PI/180,
00631 cpl_table_get_float(out, sex_mag, i, NULL),
00632 cpl_table_get_float(out, sex_magerr, i, NULL),
00633 cpl_table_get_float(out, "CLASS_STAR", i, NULL));
00634
00635 assure( !cpl_error_get_code(), return NULL,
00636 "Could not read SExtractor output table %s",
00637 filename_cat);
00638
00639 fors_star_list_insert(stars, s);
00640 }
00641 else {
00642 cpl_msg_debug(cpl_func,
00643 "Rejecting source with background = %f ADU",
00644 cpl_table_get_float(out, "BACKGROUND", i, NULL));
00645 bkg_rejected += 1;
00646 rejected += 1;
00647 }
00648 } else {
00649 rejected += 1;
00650 }
00651 }
00652
00653 cpl_msg_info(cpl_func, "%d sources sextracted, %d rejected",
00654 fors_star_list_size(stars) + rejected,
00655 rejected);
00656 }
00657
00658 if (extracted_sources != NULL) {
00659 *extracted_sources = cpl_table_duplicate(out);
00660 }
00661
00662 cleanup;
00663 return stars;
00664 }
00665
00666 #undef cleanup
00667 #define cleanup
00668
00681 static fors_star_list *
00682 extract_test(fors_extract_sky_stats *sky_stats,
00683 cpl_image **background,
00684 cpl_table **extracted_sources)
00685 {
00686 assure( sky_stats != NULL, return NULL, NULL );
00687 assure( background != NULL, return NULL, NULL );
00688
00689 sky_stats->mean = 1;
00690 sky_stats->median = 2;
00691 sky_stats->rms = 3;
00692
00693 *background = cpl_image_new(10, 20, CPL_TYPE_FLOAT);
00694
00695 fors_star_list *stars = fors_star_list_new();
00696
00697 struct {
00698 double x, y, magnitude, dmagnitude;
00699 }
00700 data[] = {
00701 {100 , 200, -10, 0.01},
00702 {1100, 200, -11, 0.01},
00703 {1 , 5 , -10, 0.01},
00704 {100 , 1200, -12, 0.01},
00705 {1100, 1200, -13, 0.01}
00706 };
00707
00708 int N = sizeof(data) / sizeof(*data);
00709 int i;
00710 double a = 2;
00711 double b = 1;
00712 double fwhm = 1.5;
00713 double orientation = 1.0;
00714
00715 for (i = 0; i < N; i++) {
00716 fors_star_list_insert(stars,
00717 fors_star_new(data[i].x,
00718 data[i].y,
00719 fwhm,
00720 a, b, orientation,
00721 data[i].magnitude,
00722 data[i].dmagnitude,
00723 1.0));
00724 }
00725
00726 if (extracted_sources != NULL) {
00727 *extracted_sources = fors_create_sources_table(stars);
00728
00729 assure (!cpl_error_get_code(), return NULL,
00730 "Could not create extracted sources table");
00731 }
00732
00733 return stars;
00734 }
00735