fors_identify.c

00001 /* $Id: fors_identify.c,v 1.40 2008/07/17 11:40:09 cizzo Exp $
00002  *
00003  * This file is part of the FORS Library
00004  * Copyright (C) 2002-2006 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
00019  */
00020 
00021 /*
00022  * $Author: cizzo $
00023  * $Date: 2008/07/17 11:40:09 $
00024  * $Revision: 1.40 $
00025  * $Name:  $
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     /* Cast is safe, see caller */
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;  /* Subset of catalog stars
00350                                             which are inside the CCD */
00351 
00352     fors_std_star_list *std_ccd_bright = NULL; /* Brightest std stars */
00353 
00354     fors_star_list *source_bright  = NULL;     /* Subset of brightest stars */
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       fors_star_print_list(CPL_MSG_ERROR, stars);
00364       fors_std_star_print_list(CPL_MSG_ERROR, cat); 
00365     */
00366     
00367 
00368     cpl_msg_info(cpl_func, "Pattern matching");
00369     cpl_msg_indent_more();
00370     
00371     /* Select standards inside CCD */
00372     {
00373         double tolerance = 100; /* pixels */
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; /* Anything can go here, not used */
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, &region);
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                 /* Searching for the nearest neighbour will permute the list,
00413                    do not do that while iterating the same list */
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                     /* If very close, remove the one with the largest magnitude
00443                        error. 
00444                        
00445                        Do not try to combine the two magnitudes because
00446                        those estimates may or may not be correlated 
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     /* Select brightest std */
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         //fors_std_star_print_list(std_ccd);
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         /* Select sources */
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     /* Finally, make the identification if a unique source is found 
00559        within the corrected position search radius.
00560 
00561        Use all catalog stars (inside CCD).
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     } /* while no identifications made */
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     /* Sketch to recompute shift:
00669      *    Get median shifts (in x and y) of identified stars.
00670      *    Then do the equivalent of a single tour through the while() loop above
00671      *    
00672      *    This could perhaps be unified with the code above to avoid duplication
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     /* Build patterns */
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; /* minimum distance between two 
00794                                     points in a pattern */
00795     double sigma_std = 0.0; /* No uncertainty of catalog patterns */
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; /* not used in this case */
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     /* Identify, 
00823        get average scale, orientation */
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         /* Make identification if patterns match within error bars
00846            (defined above as sigma_source) 
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     /* no source patterns to match */
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         /* Represent each angle as a unit vector, compute average unit vector orientation.
00875            Compute median absolute deviation with respect to this average 
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             /* To model any orientation, we should use higher order (than zero) polynomials
00895                for the shifts */
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         /* Average distance to nearest false match */
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         /* Get shift from matches */
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                     /* Compute shift of the two patterns' reference stars */
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         /* If this happens it is likely a bug, because
00956            we already checked that 'scales' is non-empty.
00957 
00958            But do not try to get the median of an empty list in any case,
00959            which would cause a crash
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     

Generated on Wed Sep 10 07:31:51 2008 for FORS Pipeline Reference Manual by  doxygen 1.4.6