/* * match: a package to match lists of stars (or other items) * Copyright (C) 2000 Michael William Richmond * * Contact: Michael William Richmond * Physics Department * Rochester Institute of Technology * 85 Lomb Memorial Drive * Rochester, NY 14623-5603 * E-mail: mwrsps@rit.edu * * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * */ /* * * FILE: match.c * * * This file contains the 'main' routine, which takes instructions * from the user, reads information from files, calls the matching * routines, and writes information back onto disk. * * Usage: match starfile1 x1 y1 mag1 starfile2 x2 y2 mag2 * [outfile=] [trirad=] [nobj=] [matchrad=] [scale=] * [max_iter=] [halt_sigma=] * [transonly] [recalc] [--version] [--help] [help] * * * * * 7/18/96 - Added transonly option. MWR * * 6/1/2000 - Added 'recalc' option; this goes through the usual steps: * * 1. read in complete lists of stars from two sources * 2. using only N brightest, find a set of matching pairs * 3. calc TRANS * 4. apply TRANS to all of list A * 5. match all of list A (with TRANS'ed coords) against * all of list B to find _many_ pairs * * but now it also adds the following * * 6. using these _many_ pairs, calc TRANS again * * The point is to avoid running the similar-triangles code * on the _many_ pairs, since that's the slow part; * we can very quickly calc a new, more accurate TRANS from * the _many_ pairs. * * MWR * * 6/14/2000 - Added "--version" and "--help" options, plus GPL. * MWR * * 7/30/2000 - Added "id1" and "id2" options. MWR * * 1/21/2001 - Added "max_iter" and "halt_sigma" options, plus some sanity * checks on values of user-supplied arguments. * MWR * * * * */ #include #include #include #include #include "misc.h" #include "match.h" #include "atpmatch.h" #undef DEBUG /* get some of diagnostic output */ static void reset_copy_ids(int numA, struct s_star *star_list_A_copy); static int reset_A_coords(int numA, struct s_star *post_list_A, struct s_star *pre_list_A); #define USAGE "usage: match starfile1 x1 y1 mag1 starfile2 x2 y2 mag2 "\ " [id1=] [id2=] [max_iter=] [halt_sigma=] "\ " [outfile=] [trirad=] [nobj=] [matchrad=] [scale=] "\ " [transonly] [recalc] [linear|quadratic|cubic] "\ " [--version] [--help] [help] " char *progname = "match"; int main ( int argc, char *argv[] ) { int i, ret; int numA, numB; int num_matched_A, num_matched_B; int numA_copy; int x1col, y1col, mag1col; int x2col, y2col, mag2col; int id1col = -1, id2col = -1; int max_iter = AT_MATCH_MAXITER; int trans_order = AT_TRANS_LINEAR; char *fileA, *fileB; double triangle_radius = AT_TRIANGLE_RADIUS; /* in triangle-space coords */ double match_radius = AT_MATCH_RADIUS; /* in units of list B */ double scale = -1.0; double halt_sigma = AT_MATCH_HALTSIGMA; int nobj = AT_MATCH_NBRIGHT; int transonly = 0; /* if 1, only find TRANS */ int recalc = 0; /* if 1, calc TRANS again */ char outfile[100]; char matched_file[100]; struct s_star *star_list_A, *star_list_B; struct s_star *star_list_A_copy; struct s_star *matched_list_A, *matched_list_B; TRANS trans; strcpy(outfile, AT_MATCH_OUTFILE); /* make special search for --version, --help, and help */ if ((argc < 2) || (strncmp(argv[1], "--help", 6) == 0) || (strncmp(argv[1], "help", 4) == 0)) { printf("%s\n", USAGE); exit(0); } if (strncmp(argv[1], "--version", 9) == 0) { printf("%s: version %s\n", progname, VERSION); exit(0); } /* * If we get here, the user seems to want to run the program normally. * Make sure there are at least the required arguments */ if (argc < 8) { fprintf(stderr, "%s\n", USAGE); exit(1); } /* these are all required args for normal usage */ fileA = argv[1]; if (sscanf(argv[2], "%d", &x1col) != 1) { shFatal("invalid argument for column for X values in first file"); } if (sscanf(argv[3], "%d", &y1col) != 1) { shFatal("invalid argument for column for Y values in first file"); } if (sscanf(argv[4], "%d", &mag1col) != 1) { shFatal("invalid argument for column for mag values in first file"); } fileB = argv[5]; if (sscanf(argv[6], "%d", &x2col) != 1) { shFatal("invalid argument for column for X values in second file"); } if (sscanf(argv[7], "%d", &y2col) != 1) { shFatal("invalid argument for column for Y values in second file"); } if (sscanf(argv[8], "%d", &mag2col) != 1) { shFatal("invalid argument for column for mag values in second file"); } /* * check for optional arguments */ for (i = 9; i < argc; i++) { if (strncmp(argv[i], "nobj=", 5) == 0) { if (sscanf(argv[i] + 5, "%d", &nobj) != 1) { shFatal("invalid argument for nobj argument"); } if (nobj <= 0) { shFatal("invalid value %d for nobj: must be > 0", nobj); } } if (strncmp(argv[i], "trirad=", 7) == 0) { if (sscanf(argv[i] + 7, "%lf", &triangle_radius) != 1) { shFatal("invalid argument for trirad argument"); } if (triangle_radius <= 0) { shFatal("invalid value %lf for triangle_radius: must be > 0", triangle_radius); } } if (strncmp(argv[i], "matchrad=", 9) == 0) { if (sscanf(argv[i] + 9, "%lf", &match_radius) != 1) { shFatal("invalid argument for matchrad argument"); } if (match_radius <= 0) { shFatal("invalid value %lf for match_radius: must be > 0", match_radius); } } if (strncmp(argv[i], "scale=", 6) == 0) { if (sscanf(argv[i] + 6, "%lf", &scale) != 1) { shFatal("invalid argument for scale argument"); } if (scale <= 0) { shFatal("invalid value %lf for scale: must be > 0", scale); } } if (strncmp(argv[i], "outfile=", 8) == 0) { if (sscanf(argv[i] + 8, "%s", outfile) != 1) { shFatal("invalid argument for outfile argument"); } } if (strncmp(argv[i], "transonly", 9) == 0) { transonly = 1; } if (strncmp(argv[i], "recalc", 6) == 0) { recalc = 1; } if (strncmp(argv[i], "linear", 6) == 0) { trans_order = AT_TRANS_LINEAR; } if (strncmp(argv[i], "quadratic", 9) == 0) { trans_order = AT_TRANS_QUADRATIC; } if (strncmp(argv[i], "cubic", 5) == 0) { trans_order = AT_TRANS_CUBIC; } if (strncmp(argv[i], "id1", 3) == 0) { if (sscanf(argv[i] + 4, "%d", &id1col) != 1) { shFatal("invalid argument for id1col argument"); } if (id1col < 0) { shFatal("invalid value %d for id1: must be >= 0", id1col); } } if (strncmp(argv[i], "id2", 3) == 0) { if (sscanf(argv[i] + 4, "%d", &id2col) != 1) { shFatal("invalid argument for id2col argument"); } if (id2col < 0) { shFatal("invalid value %d for id2col: must be >= 0", id2col); } } if (strncmp(argv[i], "max_iter", 8) == 0) { if (sscanf(argv[i] + 9, "%d", &max_iter) != 1) { shFatal("invalid argument for max_iter argument"); } if (max_iter < 0) { shFatal("invalid value %d for max_iter: must be >= 0", max_iter); } } if (strncmp(argv[i], "halt_sigma", 10) == 0) { if (sscanf(argv[i] + 11, "%lf", &halt_sigma) != 1) { shFatal("invalid argument for halt_sigma argument"); } if (halt_sigma <= 0) { shFatal("invalid value %lf for halt_sigma: must be > 0", halt_sigma); } } } #ifdef DEBUG printf("using trans_order %d\n", trans_order); #endif /* set the order used for all TRANSformations */ atTransOrderSet(trans_order); trans.order = trans_order; /* read information from the first file */ if (read_star_file(fileA, x1col, y1col, mag1col, id1col, -1, &numA, &star_list_A) != SH_SUCCESS) { shError("can't read data from file %s", fileA); exit(1); } /* * If the user specifies the 'recalc' option, we'll need to * make a second calculation of the TRANS, using the complete * set of matched pairs. We'll use this copy of the stars in set A * to replace the output matched coords (which have been converted * to those in set B) with the original coords. */ if (recalc == 1) { if (read_star_file(fileA, x1col, y1col, mag1col, id1col, -1, &numA_copy, &star_list_A_copy) != SH_SUCCESS) { shError("can't read data from file %s", fileA); exit(1); } /* sanity check */ shAssert(numA_copy == numA); /* * reset the 'id' field values in the star_list_A_copy so that they * match the 'id' field values in their counterparts in star_list_A */ reset_copy_ids(numA, star_list_A_copy); } /* read information from the second file */ if (read_star_file(fileB, x2col, y2col, mag2col, id2col, -1, &numB, &star_list_B) != SH_SUCCESS) { shError("can't read data from file %s", fileB); exit(1); } /* call the matching routines */ ret = atFindTrans(numA, star_list_A, numB, star_list_B, triangle_radius, nobj, scale, max_iter, halt_sigma, &trans); /* * if the "transonly" flag is set, stop at this point */ if (transonly == 1) { print_trans(&trans); return(0); } /* * having found the TRANS that takes A -> B, let us apply * it to all the elements in A; thus, we'll have two sets of * of stars, each in the same coordinate system */ atApplyTrans(numA, star_list_A, &trans); /* * now match up the two sets of items, and find four subsets: * * those from list A that do have matches in list B * those from list B that do have matches in list A * those from list A that do NOT have matches in list B * those from list B that do NOT have matches in list A * */ atMatchLists(numA, star_list_A, numB, star_list_B, match_radius, outfile); /* * if the user requests the 'recalc' option, then we must re-set * the coords of stars from list A to their original values; * then we can use the matched elements of lists A and B to * define an improved TRANS */ if (recalc == 1) { /* * read in new lists of stars from the output "matched" files. * The matched files have names based on 'outfile', but with * extensions "mtA" and "mtB". Ex: if 'outfile' is "matched", * then the two lists are in * * matched.mtA matched.mtB * * The format of each file is one star per line, with 4 fields: * * internal_ID xval yval mag * * where the coords (xval, yval) are in system of list B. */ sprintf(matched_file, "%s.mtA", outfile); if (read_matched_file(matched_file, &num_matched_A, &matched_list_A) != SH_SUCCESS) { shError("read_matched_file can't read data from file %s", matched_file); exit(1); } sprintf(matched_file, "%s.mtB", outfile); if (read_matched_file(matched_file, &num_matched_B, &matched_list_B) != SH_SUCCESS) { shError("read_matched_file can't read data from file %s", matched_file); exit(1); } if (reset_A_coords(num_matched_A, matched_list_A, star_list_A_copy) != 0) { fprintf(stderr, "%s: reset_A_coords returns with error\n", progname); exit(1); } /* okay, now we're ready to call atRecalcTrans */ if (atRecalcTrans(num_matched_A, matched_list_A, num_matched_B, matched_list_B, max_iter, halt_sigma, &trans) != SH_SUCCESS) { fprintf(stderr, "%s: atRecalcTrans fails\n", progname); exit(1); } } print_trans(&trans); return(0); } /*********************************************************************** * FUNCTION: reset_copy_ids * * Modify the 'id' field values in the given list (a copy of list A) * so that they will match the 'id' values in the corresponding * stars of list A. * * We have to do this because the routine which creates new s_star * structs keeps incrementing the 'id' values, and so the copies * have a different value. * * RETURNS * nothing */ static void reset_copy_ids ( int numA, /* I: number of stars in list A */ struct s_star *star_list_A_copy /* I/O: copy of star A list */ ) { int i; struct s_star *star; for (i = 0, star = star_list_A_copy; i < numA; i++, star = star->next) { shAssert(star != NULL); star->id -= numA; } } /*********************************************************************** * FUNCTION: reset_A_coords * * Given the number of elements in list 'post_list_A', * and two versions of the * stars in list A (after conversion to coord system of list B, * and before conversion), restore the original coords of stars * in the the matched list. * * RETURNS * 0 if all goes well * 1 otherwise */ static int reset_A_coords ( int numA, /* I: number of stars in the list */ struct s_star *post_list_A, /* I/O: stars in A, after they have */ /* been matched to stars in */ /* list B; coords have been */ /* converted. We'll reset */ /* coords in this list */ struct s_star *pre_list_A /* I: stars in A, with original coords */ ) { int post_index; int found_it; struct s_star *pre_star, *post_star; /* sanity checks */ shAssert(post_list_A != NULL); shAssert(pre_list_A != NULL); for (post_index = 0, post_star = post_list_A; post_index < numA; post_index++, post_star = post_star->next) { shAssert(post_star != NULL); found_it = 0; pre_star = pre_list_A; while (pre_star != NULL) { if (pre_star->id == post_star->id) { post_star->x = pre_star->x; post_star->y = pre_star->y; found_it = 1; break; } pre_star = pre_star->next; } if (found_it == 0) { printf("reset_A_coords: no match for post_star %d?\n", post_index); shAssert(0); } } return(0); }