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_identify.h>
00033
00034 #include <fors_image.h>
00035 #include <fors_pattern.h>
00036 #include <fors_point.h>
00037 #include <fors_dfs.h>
00038 #include <fors_utils.h>
00039 #include <fors_double.h>
00040
00041 #include <cpl.h>
00042
00043 #include <math.h>
00044 #include <assert.h>
00045
00052 struct _identify_method
00053 {
00054 int ncat;
00055 double nsource;
00056 double kappa;
00057 double search;
00058 double max_search;
00059 };
00060
00061 static bool
00062 inside_region(const fors_std_star *std,
00063 void *reg);
00064
00065 static void
00066 match_patterns(const fors_star_list *stars,
00067 const fors_std_star_list *std,
00068 double kappa,
00069 double *sx_00,
00070 double *sy_00,
00071 double *med_scale,
00072 double *med_angle);
00073
00079 void
00080 fors_identify_define_parameters(cpl_parameterlist *parameters,
00081 const char *context)
00082 {
00083 cpl_parameter *p;
00084 const char *full_name = NULL;
00085 const char *name;
00086
00087 name = "ncat";
00088 full_name = cpl_sprintf("%s.%s", context, name);
00089 p = cpl_parameter_new_value(full_name,
00090 CPL_TYPE_INT,
00091 "Number of catalog stars to use in "
00092 "pattern matching",
00093 context,
00094 10);
00095
00096 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00097 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00098 cpl_parameterlist_append(parameters, p);
00099 cpl_free((void *)full_name);
00100
00101
00102 name = "nsource";
00103 full_name = cpl_sprintf("%s.%s", context, name);
00104 p = cpl_parameter_new_value(full_name,
00105 CPL_TYPE_DOUBLE,
00106 "Number of sources to use in "
00107 "pattern matching, pr. catalog star",
00108 context,
00109 3.0);
00110
00111 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00112 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00113 cpl_parameterlist_append(parameters, p);
00114 cpl_free((void *)full_name);
00115
00116
00117 name = "kappa";
00118 full_name = cpl_sprintf("%s.%s", context, name);
00119 p = cpl_parameter_new_value(full_name,
00120 CPL_TYPE_DOUBLE,
00121 "Rejection parameter (scale, angle) used "
00122 "in pattern matching ",
00123 context,
00124 5.0);
00125 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00126 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00127 cpl_parameterlist_append(parameters, p);
00128 cpl_free((void *)full_name);
00129
00130
00131 name = "search";
00132 full_name = cpl_sprintf("%s.%s", context, name);
00133 p = cpl_parameter_new_value(full_name,
00134 CPL_TYPE_DOUBLE,
00135 "Search radius (pixels)",
00136 context,
00137 20.0);
00138 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00139 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00140 cpl_parameterlist_append(parameters, p);
00141 cpl_free((void *)full_name);
00142
00143 name = "maxsearch";
00144 full_name = cpl_sprintf("%s.%s", context, name);
00145 p = cpl_parameter_new_value(full_name,
00146 CPL_TYPE_DOUBLE,
00147 "Maximum search radius (pixels)",
00148 context,
00149 20.0);
00150 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00151 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00152 cpl_parameterlist_append(parameters, p);
00153 cpl_free((void *)full_name);
00154
00155 return;
00156 }
00157
00158 #undef cleanup
00159 #define cleanup \
00160 do { \
00161 cpl_free((void *)name); \
00162 } while (0)
00163
00172 identify_method *
00173 fors_identify_method_new(const cpl_parameterlist *parameters, const char *context)
00174 {
00175 identify_method *im = cpl_malloc(sizeof(*im));
00176 const char *name = NULL;
00177
00178 cpl_msg_info(cpl_func, "Identification parameters:");
00179
00180 cpl_msg_indent_more();
00181 name = cpl_sprintf("%s.%s", context, "ncat");
00182 im->ncat = dfs_get_parameter_int_const(parameters,
00183 name);
00184 cpl_free((void *)name); name = NULL;
00185 cpl_msg_indent_less();
00186
00187
00188 cpl_msg_indent_more();
00189 name = cpl_sprintf("%s.%s", context, "nsource");
00190 im->nsource = dfs_get_parameter_double_const(parameters,
00191 name);
00192 cpl_free((void *)name); name = NULL;
00193 cpl_msg_indent_less();
00194
00195
00196 cpl_msg_indent_more();
00197 name = cpl_sprintf("%s.%s", context, "kappa");
00198 im->kappa = dfs_get_parameter_double_const(parameters,
00199 name);
00200 cpl_free((void *)name); name = NULL;
00201 cpl_msg_indent_less();
00202
00203
00204 cpl_msg_indent_more();
00205 name = cpl_sprintf("%s.%s", context, "search");
00206 im->search = dfs_get_parameter_double_const(parameters,
00207 name);
00208 cpl_free((void *)name); name = NULL;
00209 cpl_msg_indent_less();
00210
00211
00212 cpl_msg_indent_more();
00213 name = cpl_sprintf("%s.%s", context, "maxsearch");
00214 im->max_search = dfs_get_parameter_double_const(parameters,
00215 name);
00216 cpl_free((void *)name); name = NULL;
00217 cpl_msg_indent_less();
00218
00219
00220 assure( !cpl_error_get_code(), return NULL, NULL );
00221
00222 return im;
00223 }
00224
00228 void
00229 fors_identify_method_delete(identify_method **em)
00230 {
00231 if (em && *em) {
00232 cpl_free(*em); *em = NULL;
00233 }
00234 return;
00235 }
00236
00246 static bool
00247 std_brighter_than(const fors_std_star *s1,
00248 void *s2)
00249 {
00250 return fors_std_star_brighter_than(s1, s2, NULL);
00251 }
00252
00262 static bool
00263 star_brighter_than(const fors_star *s1,
00264 void *s2)
00265 {
00266 return fors_star_brighter_than(s1, s2, NULL);
00267 }
00268
00277 static double
00278 distsq_shift(const fors_star *s,
00279 const fors_std_star *std,
00280 double shiftx,
00281 double shifty)
00282 {
00283 fors_point *shifted_pos = fors_point_new(std->pixel->x + shiftx,
00284 std->pixel->y + shifty);
00285
00286 double result = fors_point_distsq(s->pixel, shifted_pos);
00287
00288 fors_point_delete(&shifted_pos);
00289
00290 return result;
00291 }
00292
00301 static bool
00302 star_nearer(const fors_star *s1,
00303 const fors_star *s2,
00304 void *data)
00305 {
00306 struct {
00307 double shift_x, shift_y;
00308 const fors_std_star *ref;
00309 } *d = data;
00310
00311
00312 return
00313 distsq_shift(s1, d->ref, d->shift_x, d->shift_y) <
00314 distsq_shift(s2, d->ref, d->shift_x, d->shift_y);
00315 }
00316
00317
00318 #undef cleanup
00319 #define cleanup \
00320 do { \
00321 fors_std_star_list_delete(&std_ccd , fors_std_star_delete); \
00322 fors_std_star_list_delete(&std_ccd_bright, fors_std_star_delete); \
00323 fors_star_list_delete(&source_bright, fors_star_delete); \
00324 } while (0)
00325
00344 void
00345 fors_identify(fors_star_list *stars,
00346 fors_std_star_list *cat,
00347 const identify_method *im)
00348 {
00349 fors_std_star_list *std_ccd = NULL;
00350
00351
00352 fors_std_star_list *std_ccd_bright = NULL;
00353
00354 fors_star_list *source_bright = NULL;
00355
00356 assure( stars != NULL, return, NULL );
00357
00358 cpl_msg_info(cpl_func, "Identifying sources");
00359 cpl_msg_indent_more();
00360
00361
00362
00363
00364
00365
00366
00367
00368 cpl_msg_info(cpl_func, "Pattern matching");
00369 cpl_msg_indent_more();
00370
00371
00372 {
00373 double tolerance = 100;
00374 struct {
00375 int xlo, ylo;
00376 int xhi, yhi;
00377 } region;
00378 if (fors_star_list_size(stars) > 0) {
00379 region.xlo = fors_star_get_x(fors_star_list_min_val(stars,
00380 fors_star_get_x,
00381 NULL), NULL) - tolerance;
00382 region.ylo = fors_star_get_y(fors_star_list_min_val(stars,
00383 fors_star_get_y,
00384 NULL), NULL) - tolerance;
00385 region.xhi = fors_star_get_x(fors_star_list_max_val(stars,
00386 fors_star_get_x,
00387 NULL), NULL) + tolerance;
00388 region.yhi = fors_star_get_y(fors_star_list_max_val(stars,
00389 fors_star_get_y,
00390 NULL), NULL) + tolerance;
00391
00392 } else {
00393 region.xlo = 1;
00394 region.ylo = 1;
00395 region.xhi = 1000;
00396 region.yhi = 1000;
00397 }
00398 cpl_msg_debug(cpl_func, "Search region = (%d, %d) - (%d, %d)",
00399 region.xlo, region.ylo,
00400 region.xhi, region.yhi);
00401
00402 std_ccd = fors_std_star_list_extract(
00403 cat, fors_std_star_duplicate, inside_region, ®ion);
00404
00405 int multiple_entries = 0;
00406
00407 if (fors_std_star_list_size(std_ccd) > 1) {
00408 bool found_double;
00409 do {
00410 found_double = false;
00411
00412
00413
00414 fors_std_star_list *tmp =
00415 fors_std_star_list_duplicate(std_ccd,
00416 fors_std_star_duplicate);
00417
00418 cpl_msg_debug(cpl_func, "%d stars left", fors_std_star_list_size(tmp));
00419
00420 fors_std_star *std;
00421
00422 for (std = fors_std_star_list_first(tmp);
00423 std != NULL && !found_double;
00424 std = fors_std_star_list_next(tmp)) {
00425
00426 fors_std_star *self = fors_std_star_list_kth_val(
00427 std_ccd, 1,
00428 (fors_std_star_list_func_eval)
00429 fors_std_star_dist_arcsec,
00430 std);
00431
00432 fors_std_star *nn = fors_std_star_list_kth_val(
00433 std_ccd, 2,
00434 (fors_std_star_list_func_eval)
00435 fors_std_star_dist_arcsec,
00436 std);
00437
00438 double min_dist = fors_std_star_dist_arcsec(std, nn);
00439
00440 cpl_msg_debug(cpl_func, "dist = %f arcseconds", min_dist);
00441
00442
00443
00444
00445
00446
00447
00448 if (min_dist < 5) {
00449 multiple_entries += 1;
00450
00451 if (std->dmagnitude > nn->dmagnitude) {
00452 fors_std_star_list_remove(std_ccd, self);
00453 fors_std_star_delete(&self);
00454 } else {
00455 fors_std_star_list_remove(std_ccd, nn);
00456 fors_std_star_delete(&nn);
00457 }
00458 found_double = true;
00459 }
00460 }
00461
00462 fors_std_star_list_delete(&tmp,
00463 fors_std_star_delete);
00464
00465 } while (found_double);
00466 }
00467
00468 cpl_msg_info(cpl_func, "Found %d catalog star%s inside detector, "
00469 "ignored %d repeated source%s",
00470 fors_std_star_list_size(std_ccd),
00471 fors_std_star_list_size(std_ccd) == 1 ? "" : "s",
00472 multiple_entries,
00473 multiple_entries == 1 ? "" : "s");
00474 }
00475
00476
00477 if (fors_std_star_list_size(std_ccd) <= im->ncat) {
00478 std_ccd_bright = fors_std_star_list_duplicate(std_ccd,
00479 fors_std_star_duplicate);
00480 }
00481 else {
00482 fors_std_star *kth =
00483 fors_std_star_list_kth(std_ccd,
00484 im->ncat + 1,
00485 fors_std_star_brighter_than, NULL);
00486
00487
00488
00489 std_ccd_bright = fors_std_star_list_extract(
00490 std_ccd, fors_std_star_duplicate,
00491 std_brighter_than, kth);
00492 }
00493
00494 double sx_00, sy_00;
00495
00496 if (fors_std_star_list_size(std_ccd_bright) < 3) {
00497
00498 cpl_msg_warning(cpl_func, "Too few catalog stars (%d) available for pattern "
00499 "matching, assuming FITS header WCS solution",
00500 fors_std_star_list_size(std_ccd_bright));
00501
00502 sx_00 = 0;
00503 sy_00 = 0;
00504 }
00505 else {
00506
00507 double med_scale, med_angle;
00508
00509 cpl_msg_info(cpl_func, "Using %d brightest standards",
00510 fors_std_star_list_size(std_ccd_bright));
00511
00512 fors_std_star_print_list(CPL_MSG_DEBUG, std_ccd_bright);
00513
00514
00515 int n_sources =
00516 (int) (fors_std_star_list_size(std_ccd_bright)*im->nsource + 0.5);
00517
00518 if (fors_star_list_size(stars) <= n_sources) {
00519 source_bright = fors_star_list_duplicate(stars,
00520 fors_star_duplicate);
00521 }
00522 else {
00523 fors_star *kth =
00524 fors_star_list_kth(stars,
00525 n_sources + 1,
00526 fors_star_brighter_than, NULL);
00527
00528 source_bright = fors_star_list_extract(
00529 stars, fors_star_duplicate,
00530 star_brighter_than, kth);
00531 }
00532
00533 cpl_msg_info(cpl_func, "Using %d brightest sources",
00534 fors_star_list_size(source_bright));
00535 fors_star_print_list(CPL_MSG_DEBUG, source_bright);
00536
00537
00538 match_patterns(source_bright, std_ccd_bright,
00539 im->kappa,
00540 &sx_00, &sy_00, &med_scale, &med_angle);
00541
00542 assure( !cpl_error_get_code(), return, "Pattern matching failed" );
00543
00544
00545 if (med_scale > 1.1 || med_scale < 0.9) {
00546 cpl_msg_warning(cpl_func, "Unexpected scale from pattern "
00547 "matching (expected 1.0): assuming FITS header WCS solution");
00548
00549 sx_00 = 0;
00550 sy_00 = 0;
00551 }
00552
00553
00554 }
00555 cpl_msg_indent_less();
00556
00557
00558
00559
00560
00561
00562
00563
00564 int number_of_ids = 0;
00565 double search_radius = im->search;
00566 bool require_unique = true;
00567
00568 while (number_of_ids == 0 && search_radius <= im->max_search) {
00569
00570 cpl_msg_info(cpl_func, "Identification");
00571
00572 cpl_msg_indent_more();
00573 cpl_msg_info(cpl_func, "Average shift = (%.2f, %.2f) pixels",
00574 sx_00,
00575 sy_00);
00576
00577 if (fabs(sx_00) > 10.0) {
00578 cpl_msg_warning(cpl_func, "Large x-shift");
00579 }
00580 if (fabs(sy_00) > 10.0) {
00581 cpl_msg_warning(cpl_func, "Large y-shift");
00582 }
00583
00584 cpl_msg_info(cpl_func, "search_radius = %.2f pixels", search_radius);
00585
00586 struct {
00587 double shift_x, shift_y;
00588 const fors_std_star *ref;
00589 } data;
00590 data.shift_x = sx_00;
00591 data.shift_y = sy_00;
00592
00593 if (fors_star_list_size(stars) > 0) {
00594 for (data.ref = fors_std_star_list_first_const(std_ccd);
00595 data.ref != NULL;
00596 data.ref = fors_std_star_list_next_const(std_ccd)) {
00597
00598 fors_star *nearest_1 = fors_star_list_kth(stars,
00599 1,
00600 star_nearer,
00601 &data);
00602
00603 fors_star *nearest_2 = fors_star_list_kth(stars,
00604 2,
00605 star_nearer,
00606 &data);
00607
00608 cpl_msg_debug(cpl_func, "Nearest candidates at distance %f and %f pixels",
00609 sqrt(distsq_shift(nearest_1, data.ref, data.shift_x, data.shift_y)),
00610 sqrt(distsq_shift(nearest_2, data.ref, data.shift_x, data.shift_y)));
00611
00612 if (distsq_shift(nearest_1, data.ref, data.shift_x, data.shift_y)
00613 <=
00614 search_radius * search_radius) {
00615
00616 if (!require_unique ||
00617 distsq_shift(nearest_2, data.ref, data.shift_x, data.shift_y)
00618 > search_radius * search_radius) {
00619
00620 cpl_msg_debug(cpl_func, "unique source inside %f pixels",
00621 search_radius);
00622
00623 nearest_1->id = fors_std_star_duplicate(data.ref);
00624 number_of_ids += 1;
00625 }
00626 }
00627 }
00628 }
00629
00630 cpl_msg_info(cpl_func, "Identified %d star%s",
00631 number_of_ids, (number_of_ids == 1) ? "" : "s");
00632
00633 if (number_of_ids == 0) {
00634
00635 if (fabs(sx_00) > 0.1 &&
00636 fabs(sy_00) > 0.1) {
00637
00638 cpl_msg_warning(cpl_func,
00639 "No identifications made, "
00640 "set shift to zero and try again");
00641 search_radius = 20;
00642 require_unique = false;
00643
00644 sx_00 = 0;
00645 sy_00 = 0;
00646 }
00647 else {
00648 require_unique = false;
00649
00650 search_radius *= 2;
00651
00652 cpl_msg_warning(cpl_func,
00653 "No identifications made, "
00654 "double search radius and try again");
00655 }
00656 }
00657
00658 cpl_msg_indent_less();
00659
00660 }
00661
00662 if (number_of_ids == 0) {
00663 cpl_msg_warning(cpl_func,
00664 "No identifications made, max search radius reached: %f pixels",
00665 im->max_search);
00666 }
00667 else {
00668
00669
00670
00671
00672
00673
00674 }
00675
00676
00677 cpl_msg_indent_less();
00678
00679 cleanup;
00680 return;
00681 }
00682
00691 static bool
00692 inside_region(const fors_std_star *std,
00693 void *reg)
00694 {
00695 const struct {
00696 int xlo, ylo;
00697 int xhi, yhi;
00698 } *region = reg;
00699
00700 return
00701 region->xlo <= std->pixel->x && std->pixel->x <= region->xhi &&
00702 region->ylo <= std->pixel->y && std->pixel->y <= region->yhi;
00703 }
00704
00705
00706
00707 #undef cleanup
00708 #define cleanup \
00709 do { \
00710 fors_point_list_delete(&std_points, fors_point_delete); \
00711 fors_point_list_delete(&source_points, fors_point_delete); \
00712 fors_pattern_list_delete(&std_patterns, fors_pattern_delete); \
00713 fors_pattern_list_delete(&source_patterns, fors_pattern_delete); \
00714 double_list_delete(&scales, double_delete); \
00715 double_list_delete(&angles, double_delete); \
00716 double_list_delete(&angle_cos, double_delete); \
00717 double_list_delete(&angle_sin, double_delete); \
00718 double_list_delete(&match_dist, double_delete); \
00719 double_list_delete(&shiftx, double_delete); \
00720 double_list_delete(&shifty, double_delete); \
00721 } while (0)
00722
00738 static void
00739 match_patterns(const fors_star_list *stars,
00740 const fors_std_star_list *std,
00741 double kappa,
00742 double *sx_00,
00743 double *sy_00,
00744 double *med_scale,
00745 double *med_angle)
00746 {
00747 fors_point_list *std_points = NULL;
00748 fors_point_list *source_points = NULL;
00749
00750 fors_pattern_list *std_patterns = NULL;
00751 fors_pattern_list *source_patterns = NULL;
00752
00753 double_list *scales = NULL;
00754 double_list *angles = NULL;
00755 double_list *angle_cos = NULL;
00756 double_list *angle_sin = NULL;
00757 double_list *match_dist = NULL;
00758 double_list *shiftx = NULL;
00759 double_list *shifty = NULL;
00760
00761 assure( sx_00 != NULL, return, NULL );
00762 assure( sy_00 != NULL, return, NULL );
00763
00764
00765 std_points = fors_point_list_new();
00766 {
00767 const fors_std_star *s;
00768
00769 for (s = fors_std_star_list_first_const(std);
00770 s != NULL;
00771 s = fors_std_star_list_next_const(std)) {
00772
00773 fors_point_list_insert(std_points,
00774 fors_point_new(s->pixel->x,
00775 s->pixel->y));
00776 }
00777 }
00778
00779 source_points = fors_point_list_new();
00780 {
00781 const fors_star *s;
00782
00783 for (s = fors_star_list_first_const(stars);
00784 s != NULL;
00785 s = fors_star_list_next_const(stars)) {
00786
00787 fors_point_list_insert(source_points,
00788 fors_point_new(s->pixel->x,
00789 s->pixel->y));
00790 }
00791 }
00792
00793 const double min_dist = 1.0;
00794
00795 double sigma_std = 0.0;
00796 std_patterns =
00797 fors_pattern_new_from_points(std_points, min_dist, sigma_std);
00798 cpl_msg_info(cpl_func, "Created %d catalog patterns",
00799 fors_pattern_list_size(std_patterns));
00800
00801 double sigma_source;
00802 if (fors_star_list_size(stars) > 0) {
00803 sigma_source = fors_star_list_median(stars,
00804 fors_star_extension, NULL);
00805 cpl_msg_info(cpl_func, "Average source extension = %.2f pixels", sigma_source);
00806 } else {
00807 sigma_source = -1;
00808 }
00809 source_patterns =
00810 fors_pattern_new_from_points(source_points, min_dist, sigma_source);
00811
00812 cpl_msg_info(cpl_func, "Created %d source patterns",
00813 fors_pattern_list_size(source_patterns));
00814
00815 scales = double_list_new();
00816 angles = double_list_new();
00817 angle_cos = double_list_new();
00818 angle_sin = double_list_new();
00819 match_dist = double_list_new();
00820
00821 if ( fors_pattern_list_size(source_patterns) > 0) {
00822
00823
00824 fors_pattern *p;
00825
00826 for (p = fors_pattern_list_first(std_patterns);
00827 p != NULL;
00828 p = fors_pattern_list_next(std_patterns)) {
00829
00830 fors_pattern *nearest_source =
00831 fors_pattern_list_min_val(source_patterns,
00832 (fors_pattern_list_func_eval) fors_pattern_distsq,
00833 p);
00834
00835 double scale = fors_pattern_get_scale(p, nearest_source);
00836 double angle = fors_pattern_get_angle(p, nearest_source);
00837 double angle_c = cos(angle);
00838 double angle_s = sin(angle);
00839 double dist = sqrt(fors_pattern_distsq(p, nearest_source));
00840 double dist_per_error = fors_pattern_dist_per_error(p, nearest_source);
00841
00842 cpl_msg_debug(cpl_func, "dist, ndist, scale, orientation = %f, %f, %f, %f",
00843 dist, dist_per_error, scale, angle * 360/(2*M_PI));
00844
00845
00846
00847
00848 if (dist_per_error < 1.0) {
00849 double_list_insert(scales , double_duplicate(&scale));
00850 double_list_insert(angles , double_duplicate(&angle));
00851 double_list_insert(angle_cos, double_duplicate(&angle_c));
00852 double_list_insert(angle_sin, double_duplicate(&angle_s));
00853 double_list_insert(match_dist, double_duplicate(&dist));
00854 }
00855 }
00856 }
00857 else {
00858
00859 }
00860
00861 if ( double_list_size(scales) >= 2 ) {
00862 double scale_avg = double_list_median(scales, double_eval, NULL);
00863 double scale_stdev = double_list_mad(scales, double_eval, NULL) * STDEV_PR_MAD;
00864
00865 cpl_msg_info(cpl_func, "Median scale = %.4f (standard deviation = %.4f)",
00866 scale_avg, scale_stdev);
00867
00868 *med_scale = scale_avg;
00869
00870 if (scale_stdev > 0.2) {
00871 cpl_msg_warning(cpl_func, "Uncertain scale determination");
00872 }
00873
00874
00875
00876
00877 double angle_avg = atan2(double_list_mean(angle_sin, double_eval, NULL),
00878 double_list_mean(angle_cos, double_eval, NULL));
00879 double angle_stdev = STDEV_PR_MAD *
00880 double_list_median(angles, (double_list_func_eval) fors_angle_diff, &angle_avg);
00881
00882 cpl_msg_info(cpl_func, "Average orientation = %.1f (standard deviation = %.1f degrees)",
00883 angle_avg * 360 / (2*M_PI),
00884 angle_stdev * 360 / (2*M_PI));
00885
00886 *med_angle = angle_avg;
00887
00888 if (angle_avg > M_PI/4 || angle_avg < -M_PI/4) {
00889 cpl_msg_warning(cpl_func, "Expected orientation = 0 degrees");
00890 }
00891
00892 if (angle_avg > M_PI/4 || angle_avg < -M_PI/4) {
00893 cpl_msg_warning(cpl_func, "Correction model assumes orientation = 0 degrees");
00894
00895
00896 }
00897
00898 double avg_dist = double_list_mean(match_dist, double_eval, NULL);
00899 double false_dist = 1.0/sqrt(M_PI * fors_pattern_list_size(source_patterns));
00900
00901
00902 cpl_msg_info(cpl_func, "Average match distance = %f", avg_dist);
00903 cpl_msg_info(cpl_func, "False match distance = %f", false_dist);
00904 cpl_msg_info(cpl_func, "Safety index = %.3f (should be >~ 5)", false_dist / avg_dist);
00905 if (false_dist / avg_dist < 1.5) {
00906 cpl_msg_warning(cpl_func, "Uncertain pattern matching");
00907 }
00908
00909
00910 shiftx = double_list_new();
00911 shifty = double_list_new();
00912 {
00913 fors_pattern *p;
00914
00915 for (p = fors_pattern_list_first(std_patterns);
00916 p != NULL;
00917 p = fors_pattern_list_next(std_patterns)) {
00918
00919 fors_pattern *nearest_source =
00920 fors_pattern_list_min_val(
00921 source_patterns,
00922 (fors_pattern_list_func_eval) fors_pattern_distsq,
00923 p);
00924
00925 double dist = sqrt(fors_pattern_distsq(p, nearest_source));
00926 double dist_per_error = fors_pattern_dist_per_error(p, nearest_source);
00927 double scale = fors_pattern_get_scale(p, nearest_source);
00928 double angle = fors_pattern_get_angle(p, nearest_source);
00929
00930 cpl_msg_debug(cpl_func, "scale, orientation, distance, norm.distance "
00931 "= %f, %f, %f, %f",
00932 scale, angle * 360/(2*M_PI), dist, dist_per_error);
00933
00934 if (dist_per_error < 1.0 &&
00935 fabs(scale - scale_avg) <= kappa * scale_stdev &&
00936 fors_angle_diff(&angle, &angle_avg) <= kappa * angle_stdev) {
00937
00938
00939 double shift_x = fors_pattern_get_ref(nearest_source)->x - fors_pattern_get_ref(p)->x;
00940 double shift_y = fors_pattern_get_ref(nearest_source)->y - fors_pattern_get_ref(p)->y;
00941
00942 cpl_msg_debug(cpl_func, "Accepted, shift = (%f, %f) pixels",
00943 shift_x, shift_y);
00944
00945 double_list_insert(shiftx, double_duplicate(&shift_x));
00946 double_list_insert(shifty, double_duplicate(&shift_y));
00947 }
00948 }
00949 }
00950
00951 if (double_list_size(shiftx) > 0) {
00952 *sx_00 = double_list_median(shiftx, double_eval, NULL);
00953 }
00954 else {
00955
00956
00957
00958
00959
00960
00961 cpl_msg_warning(cpl_func, "No star identifications!");
00962 }
00963
00964 if (double_list_size(shiftx) > 0) {
00965 *sy_00 = double_list_median(shifty, double_eval, NULL);
00966 }
00967 else {
00968 cpl_msg_warning(cpl_func, "No star identifications!");
00969 }
00970 }
00971 else {
00972 cpl_msg_warning(cpl_func,
00973 "Too few (%d) matching patterns: assuming zero shift",
00974 double_list_size(scales));
00975 *sx_00 = 0;
00976 *sy_00 = 0;
00977 *med_scale = 1.0;
00978 *med_angle = 0.0;
00979 }
00980 cpl_msg_indent_less();
00981
00982 cleanup;
00983 return;
00984 }
00985
00986
00987
00988