1
2 /*
3 * match: a package to match lists of stars (or other items)
4 * Copyright (C) 2000 Michael William Richmond
5 *
6 * Contact: Michael William Richmond
7 * Physics Department
8 * Rochester Institute of Technology
9 * 85 Lomb Memorial Drive
10 * Rochester, NY 14623-5603
11 * E-mail: mwrsps@rit.edu
12 *
13 *
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
18 *
19 * This program is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with this program; if not, write to the Free Software
26 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
27 *
28 */
29
30 /*
31 * <AUTO>
32 * FILE: atpmatch.c
33 *
34 * <HTML>
35 * This file contains routines that try to match up items
36 * in two different lists, which might have very different
37 * coordinate systems.
38 *
39 * Stars must have been placed into "s_star" structures before
40 * being passed to the functions in this file. Note that
41 * the "x" and "y" fields of an s_star may contain (RA, Dec),
42 * or (row, col), or some other coordinates.
43 * </HTML>
44 *
45 * </AUTO>
46 *
47 */
48
49 /*
50 * -------------------------------------------------------------------------
51 * atFindTrans public Find TRANS to match coord systems of
52 * two lists of items
53 * atApplyTrans public Apply a TRANS to a list of items,
54 * modifying two of the elements in each item
55 * atMatchLists public Find all pairs of items on two lists
56 * which match (and don't match)
57 * atRecalcTrans public Calculate a TRANS, given two lists of
58 * stars which _already_ have been matched
59 * atFindMedtf public Calculates MEDTF statistics, assuming
60 * that two lists have same scale, rotation
61 *
62 * All public functions appear at the start of this source-code file.
63 * "Private" static functions appear following them.
64 *
65 * Conditional compilation is controlled by the following macros:
66 *
67 * DEBUG
68 * DEBUG2
69 *
70 * AUTHORS: SHIVA Creation date: Jan 22, 1996
71 * Michael Richmond
72 *
73 * Modified for stand-alone Linux operation: April 26, 1996
74 * Michael Richmond
75 *
76 * Added more digits to printout in "print_trans": Aug 18, 1996
77 * Michael Richmond
78 *
79 * Changed 'iter_trans()' to reject 3-sigma outliers,
80 * instead of 2-sigma outliers. Yields many more matches
81 * at end of iterating, and works better for TASS Mark IV
82 * data. May 25, 2000
83 * Michael Richmond
84 *
85 * Added new public function "atRecalcTrans", which uses
86 * existing matched lists of stars to calculate a TRANS.
87 * Allows us to use _all_ the stars in both lists,
88 * not just the brightest N. June 1, 2000
89 * Michael Richmond
90 *
91 * Added 'recalc_flag' argument to 'iter_trans()', to allow
92 * different behavior on first iteration if we call it
93 * from atRecalcTrans or not. June 2, 2000
94 * Michael Richmond
95 *
96 * Changed the "3" in "3-sigma" outliers rejected in iter_trans
97 * into a #define value AT_MATCH_NSIGMA, defined in the
98 * .h file. June 2, 2000
99 * Michael Richmond
100 *
101 * Modified so that this single file contains routines to handle
102 * the linear, quadratic, and cubic transformation cases.
103 * June 10, 2000
104 * Michael Richmond
105 *
106 * Replaced old 'gauss_jordon' routine to solve matrix equation
107 * with new 'gauss_matrix' routine; uses Gaussian elimination
108 * with back-substitution instead of Gauss-Jordon technique.
109 * Not associated with "Numerical Recipes" in any way -- hah!
110 * June 19, 2000
111 * Michael Richmond
112 *
113 * Added MEDTF calculations and TRANS diagnostics, as suggested
114 * by John Blakeslee.
115 * Dec 12, 2001
116 * Michael Richmond
117 *
118 * Fixed off-by-one bug in "remove_elem", as suggested by
119 * Andrew Bennett.
120 * Also fixed small error in location of paranthesis in
121 * the "gauss_matrix" routine, again thanks to Andrew.
122 * Dec 28, 2001
123 * Michael Richmond
124 *
125 * Fixed bug in "atCalcRMS" which caused assertion to fail
126 * when the routine was given two empty lists.
127 * Similar bugs addressed by making early checks to the
128 * number of stars in the passed lists in other funcs:
129 * atFindMedtf
130 * Nov 22, 2002
131 * Michael Richmond
132 *
133 * Added the ratio of triangle sizes to a debugging printf
134 * statement in make_vote_matrix.
135 * June 26, 2010
136 * Michael Richmond
137 */
138
139
140 #include <stdio.h>
141 #include <math.h> /* need this for 'sqrt' in calc_distances */
142 #include <stdlib.h>
143 #include <string.h>
144 #include <math.h>
145 #include "misc.h"
146 #include "atpmatch.h"
147
148 #undef DEBUG /* get some of diagnostic output */
149 #undef DEBUG2 /* get LOTS more diagnostic output */
150 #undef DEBUG3 /* run 'test_routine' to test matrix inversion */
151
152
153 /*
154 * used in the 'gauss_matrix' matrix inversion routine,
155 * as a check for very very small numbers which might cause
156 * the matrix solution to be unstable.
157 */
158 #define MATRIX_TOL 1.0e-12
159
160 /*
161 * To evaluate the quality of a match between two sets of stars,
162 * we look at the differences in their positions after transforming
163 * those in list A to the coordinate system of list B. We sort
164 * those distances and pick the one closest to this percentile
165 * to characterize the distribution. One stdev should include
166 * about 68% of the data.
167 * This is used in routine 'iter_trans'.
168 */
169 #define ONE_STDEV_PERCENTILE 0.683
170
171 /*
172 * these values are used to tell iter_trans() whether it is being
173 * called from atRecalcTrans or not. If yes, then it is safe
174 * to include _all_ the matched pairs in finding TRANS, in the
175 * very first iteration. If no, then only use the best AT_MATCH_STARTN
176 * pairs in the first iteration.
177 */
178 #define RECALC_YES 1
179 #define RECALC_NO 0
180
181
182 /*
183 * the following are "private" functions, used internally only.
184 */
185
186 /* this typedef is used several sorting routines */
187 typedef int (*PFI)();
188
189 static int set_star(s_star *star, double x, double y, double mag);
190 static void copy_star(s_star *from_ptr, s_star *to_ptr);
191 static void copy_star_array(s_star *from_array, s_star *to_array, int num);
192 static void free_star_array(s_star *array);
193 #ifdef DEBUG
194 static void print_star_array(s_star *array, int num);
195 #endif
196 static double **calc_distances(s_star *star_array, int numstars);
197 static void free_distances(double **array, int num);
198 #ifdef DEBUG
199 static void print_dist_matrix(double **matrix, int num);
200 #endif
201 static void set_triangle(s_triangle *triangle, s_star *star_array,
202 int i, int j, int k, double **dist_matrix);
203 #ifdef DEBUG2
204 static void print_triangle_array(s_triangle *t_array, int numtriangles,
205 s_star *star_array, int numstars);
206 #endif
207 static s_triangle *stars_to_triangles(s_star *star_array, int numstars,
208 int nbright, int *numtriangles);
209 static void sort_star_by_mag(s_star *array, int num);
210 static int compare_star_by_mag(s_star *star1, s_star *star2);
211 static void sort_star_by_x(s_star *array, int num);
212 static int compare_star_by_x(s_star *star1, s_star *star2);
213 static void sort_star_by_match_id(s_star *array, int num);
214 static int compare_star_by_match_id(s_star *star1, s_star *star2);
215 static int fill_triangle_array(s_star *star_array, int numstars,
216 double **dist_matrix,
217 int numtriangles, s_triangle *t_array);
218 static void sort_triangle_array(s_triangle *array, int num);
219 static int compare_triangle(s_triangle *triangle1, s_triangle *triangle2);
220 static int find_ba_triangle(s_triangle *array, int num, double ba0);
221 static void prune_triangle_array(s_triangle *t_array, int *numtriangles);
222 static int **make_vote_matrix(s_star *star_array_A, int num_stars_A,
223 s_star *star_array_B, int num_stars_B,
224 s_triangle *t_array_A, int num_triangles_A,
225 s_triangle *t_array_B, int num_triangles_B,
226 int nbright, double radius,
227 double min_scale, double max_scale,
228 double rotation_deg, double tolerance_deg);
229 #ifdef DEBUG
230 static void print_vote_matrix(int **vote_matrix, int numcells);
231 #endif
232 static int top_vote_getters(int **vote_matrix, int num, int **winner_votes,
233 int **winner_index_A, int **winner_index_B);
234 static int calc_trans(int nbright, s_star *star_array_A, int num_stars_A,
235 s_star *star_array_B, int num_stars_B,
236 int *winner_votes,
237 int *winner_index_A, int *winner_index_B,
238 TRANS *trans);
239 static s_star *list_to_array(int num_stars, struct s_star *list);
240 static void reset_array_ids(struct s_star *star_list, int num_stars,
241 struct s_star *star_array);
242
243
244 /*
245 * these are functions used to solve a matrix equation which
246 * gives us the transformation from one coord system to the other
247 */
248 static int gauss_matrix(double **matrix, int num, double *vector);
249 static int gauss_pivot (double **matrix, int num, double *vector,
250 double *biggest_val, int row);
251
252 static double ** alloc_matrix(int n);
253 static void free_matrix(double **matrix, int n);
254 #ifdef DEBUG
255 static void print_matrix(double **matrix, int n);
256 #endif
257
258 #ifdef DEBUG3
259 static void test_routine (void);
260 #endif
261
262 static int iter_trans(int nbright, s_star *star_array_A, int num_stars_A,
263 s_star *star_array_B, int num_stars_B,
264 int *winner_votes,
265 int *winner_index_A, int *winner_index_B,
266 int recalc_flag,
267 int max_iter, double halt_sigma,
268 TRANS *trans);
269 static int compare_double(double *f1, double *f2);
270 static double find_percentile(double *array, int num, double perc);
271 static int calc_trans_coords(s_star *star, TRANS *trans,
272 double *newx, double *newy);
273 static int apply_trans(s_star *star_array, int num_stars, TRANS *trans);
274 static int double_sort_by_match_id(s_star *star_array_A, int num_stars_A,
275 s_star *star_array_B, int num_stars_B);
276 static int match_arrays_slow(s_star *star_array_A, int num_stars_A,
277 s_star *star_array_B, int num_stars_B,
278 double radius,
279 s_star **star_array_J, int *num_stars_J,
280 s_star **star_array_K, int *num_stars_K,
281 s_star **star_array_L, int *num_stars_L,
282 s_star **star_array_M, int *num_stars_M);
283 static int add_element(s_star *new_star, s_star **star_array, int *total_num,
284 int *current_num);
285 static void remove_elem(s_star *star_array, int num, int *num_stars);
286 static int remove_repeated_elements(s_star *star_array_1, int *num_stars_1,
287 s_star *star_array_2, int *num_stars_2);
288 static void remove_same_elements(s_star *star_array_1, int num_stars_1,
289 s_star *star_array_2, int *num_stars_2);
290 static void write_array(int num_stars, struct s_star *star_array,
291 char *filename);
292 static int write_small_arrays(double ra, double dec,
293 int num_stars, struct s_star *star_array,
294 int nbright, int num_triangles,
295 struct s_triangle *t_array, char *outfile);
296 static int is_desired_rotation(struct s_triangle *tri_1,
297 struct s_triangle *tri_2,
298 double want_angle_deg,
299 double tolerance_deg,
300 double *actual_angle_deg);
301
302 /*
303 * we have three different versions of a routine to do the
304 * dirty work of calculating the TRANS which best turns
305 * coords of system A into coords of system B.
306 * There is a version for linear transformations, quadratic ones,
307 * and cubic ones.
308 *
309 * All are called by the 'calc_trans' function; it uses the
310 * trans->order value to figure out which one is appropriate.
311 */
312 static int calc_trans_linear(int nbright,
313 s_star *star_array_A, int num_stars_A,
314 s_star *star_array_B, int num_stars_B,
315 int *winner_votes,
316 int *winner_index_A, int *winner_index_B,
317 TRANS *trans);
318 static int calc_trans_quadratic(int nbright,
319 s_star *star_array_A, int num_stars_A,
320 s_star *star_array_B, int num_stars_B,
321 int *winner_votes,
322 int *winner_index_A, int *winner_index_B,
323 TRANS *trans);
324 static int calc_trans_cubic(int nbright,
325 s_star *star_array_A, int num_stars_A,
326 s_star *star_array_B, int num_stars_B,
327 int *winner_votes,
328 int *winner_index_A, int *winner_index_B,
329 TRANS *trans);
330
331
332
333 /************************************************************************
334 * <AUTO EXTRACT>
335 *
336 * ROUTINE: atFindTrans
337 *
338 * DESCRIPTION:
339 * This function is based on the algorithm described in Valdes et al.,
340 * PASP 107, 1119 (1995). It tries to
341 * a. match up objects in the two chains
342 * a. find a coordinate transformation that takes coords in
343 * objects in chainA and changes to those in chainB.
344 *
345 *
346 * Actually, this is a top-level "driver" routine that calls smaller
347 * functions to perform actual tasks. It mostly creates the proper
348 * inputs and outputs for the smaller routines.
349 *
350 * RETURN:
351 * SH_SUCCESS if all goes well
352 * SH_GENERIC_ERROR if an error occurs
353 *
354 * </AUTO>
355 */
356
357 int
atFindTrans(int numA,struct s_star * listA,int numB,struct s_star * listB,double radius,int nobj,double min_scale,double max_scale,double rotation_deg,double tolerance_deg,int max_iter,double halt_sigma,TRANS * trans)358 atFindTrans
359 (
360 int numA, /* I: number of stars in list A */
361 struct s_star *listA, /* I: match this set of objects with list B */
362 int numB, /* I: number of stars in list B */
363 struct s_star *listB, /* I: match this set of objects with list A */
364 double radius, /* I: max radius in triangle-space allowed for */
365 /* a pair of triangles to match */
366 int nobj, /* I: max number of bright stars to use in creating */
367 /* triangles for matching from each list */
368 double min_scale, /* I: minimum permitted relative scale factor */
369 /* if -1, any scale factor is allowed */
370 double max_scale, /* I: maximum permitted relative scale factor */
371 /* if -1, any scale factor is allowed */
372 double rotation_deg, /* I: desired relative angle of coord systems (deg) */
373 /* if AT_MATCH_NOANGLE, any orientation is allowed */
374 double tolerance_deg, /* I: allowed range of orientation angles (deg) */
375 /* if AT_MATCH_NOANGLE, any orientation is allowed */
376 int max_iter, /* I: go through at most this many iterations */
377 /* in the iter_trans() loop. */
378 double halt_sigma, /* I: halt the fitting procedure if the mean */
379 /* residual becomes this small */
380 TRANS *trans /* O: place into this TRANS structure's fields */
381 /* the coeffs which convert coords of chainA */
382 /* into coords of chainB system. */
383 )
384 {
385 int i, nbright, min;
386 int num_stars_A; /* number of stars in chain A */
387 int num_stars_B; /* number of stars in chain B */
388 int num_triangles_A; /* number of triangles formed from chain A */
389 int num_triangles_B; /* number of triangles formed from chain B */
390 int **vote_matrix;
391 int *winner_votes; /* # votes gotten by top pairs of matched stars */
392 int *winner_index_A; /* elem i in this array is index in star array A */
393 /* which matches ... */
394 int *winner_index_B; /* elem i in this array, index in star array B */
395 int start_pairs;
396 s_star *star_array_A = NULL;
397 s_star *star_array_B = NULL;
398 s_triangle *triangle_array_A = NULL;
399 s_triangle *triangle_array_B = NULL;
400
401 num_stars_A = numA;
402 num_stars_B = numB;
403 star_array_A = list_to_array(numA, listA);
404 star_array_B = list_to_array(numB, listB);
405
406 #ifdef DEBUG3
407 test_routine();
408 #endif
409
410 shAssert(star_array_A != NULL);
411 shAssert(star_array_B != NULL);
412
413 switch (trans->order) {
414 case AT_TRANS_LINEAR:
415 start_pairs = AT_MATCH_STARTN_LINEAR;
416 break;
417 case AT_TRANS_QUADRATIC:
418 start_pairs = AT_MATCH_STARTN_QUADRATIC;
419 break;
420 case AT_TRANS_CUBIC:
421 start_pairs = AT_MATCH_STARTN_CUBIC;
422 break;
423 default:
424 shError("atFindTrans: invalid trans->order %d ", trans->order);
425 break;
426 }
427
428
429
430 /*
431 * here we check to see if each list of stars contains a
432 * required minimum number of stars. If not, we return with
433 * an error message, and SH_GENERIC_ERROR.
434 *
435 * In addition, we check to see that each list has at least 'nobj'
436 * items. If not, we set 'nbright' to the minimum of the two
437 * list lengths, and print a warning message so the user knows
438 * that we're using fewer stars than he asked.
439 *
440 * On the other hand, if the user specifies a value of "nobj" which
441 * is too SMALL, then we ignore it and use the smallest valid
442 * value (which is start_pairs).
443 */
444 min = (num_stars_A < num_stars_B ? num_stars_A : num_stars_B);
445 if (min < start_pairs) {
446 shError("atFindTrans: only %d stars in list(s), require at least %d",
447 min, start_pairs);
448 free_star_array(star_array_A);
449 free_star_array(star_array_B);
450 return(SH_GENERIC_ERROR);
451 }
452 if (nobj > min) {
453 shDebug(AT_MATCH_ERRLEVEL,
454 "atFindTrans: using only %d stars, fewer than requested %d",
455 min, nobj);
456 nbright = min;
457 }
458 else {
459 nbright = nobj;
460 }
461 if (nbright < start_pairs) {
462 shDebug(AT_MATCH_ERRLEVEL,
463 "atFindTrans: must use %d stars, more than requested %d",
464 start_pairs, nobj);
465 nbright = start_pairs;
466 }
467
468 /* this is a sanity check on the above checks */
469 shAssert((nbright >= start_pairs) && (nbright <= min));
470
471
472 #ifdef DEBUG
473 printf("here comes star array A\n");
474 print_star_array(star_array_A, num_stars_A);
475 printf("here comes star array B\n");
476 print_star_array(star_array_B, num_stars_B);
477 #endif
478
479 /*
480 * we now convert each list of stars into a list of triangles,
481 * using only a subset of the "nbright" brightest items in each list.
482 */
483 triangle_array_A = stars_to_triangles(star_array_A, num_stars_A,
484 nbright, &num_triangles_A);
485 shAssert(triangle_array_A != NULL);
486 triangle_array_B = stars_to_triangles(star_array_B, num_stars_B,
487 nbright, &num_triangles_B);
488 shAssert(triangle_array_B != NULL);
489
490
491 /*
492 * Now we prune the triangle arrays to eliminate those with
493 * ratios (b/a) > AT_MATCH_RATIO,
494 * since Valdes et al. say that this speeds things up and eliminates
495 * lots of closely-packed triangles.
496 */
497 prune_triangle_array(triangle_array_A, &num_triangles_A);
498 prune_triangle_array(triangle_array_B, &num_triangles_B);
499 #ifdef DEBUG2
500 printf("after pruning, here comes triangle array A\n");
501 print_triangle_array(triangle_array_A, num_triangles_A,
502 star_array_A, num_stars_A);
503 printf("after pruning, here comes triangle array B\n");
504 print_triangle_array(triangle_array_B, num_triangles_B,
505 star_array_B, num_stars_B);
506 #endif
507
508
509 /*
510 * Next, we want to try to match triangles in the two arrays.
511 * What we do is to create a "vote matrix", which is a 2-D array
512 * with "nbright"-by-"nbright" cells. The cell with
513 * coords [i][j] holds the number of matched triangles in which
514 *
515 * item [i] in star_array_A matches item [j] in star_array_B
516 *
517 * We'll use this "vote_matrix" to figure out a first guess
518 * at the transformation between coord systems.
519 *
520 * Note that if there are fewer than "nbright" stars
521 * in either list, we'll still make the vote_matrix
522 * contain "nbright"-by-"nbright" cells ...
523 * there will just be a lot of cells filled with zero.
524 */
525 vote_matrix = make_vote_matrix(star_array_A, num_stars_A,
526 star_array_B, num_stars_B,
527 triangle_array_A, num_triangles_A,
528 triangle_array_B, num_triangles_B,
529 nbright, radius, min_scale, max_scale,
530 rotation_deg, tolerance_deg);
531
532
533 /*
534 * having made the vote_matrix, we next need to pick the
535 * top 'nbright' vote-getters. We call 'top_vote_getters'
536 * and are given, in its output arguments, pointers to three
537 * arrays, each of which has 'nbright' elements pertaining
538 * to a matched pair of STARS:
539 *
540 * winner_votes[] number of votes of winners, in descending order
541 * winner_index_A[] index of star in star_array_A
542 * winner_index_B[] index of star in star_array_B
543 *
544 * Thus, the pair of stars which matched in the largest number
545 * of triangles will be
546 *
547 * star_array_A[winner_index_A[0]] from array A
548 * star_array_B[winner_index_A[0]] from array B
549 *
550 * and the pair of stars which matched in the second-largest number
551 * of triangles will be
552 *
553 * star_array_A[winner_index_A[1]] from array A
554 * star_array_B[winner_index_A[1]] from array B
555 *
556 * and so on.
557 */
558 top_vote_getters(vote_matrix, nbright,
559 &winner_votes, &winner_index_A, &winner_index_B);
560
561 /*
562 * here, we disqualify any of the top vote-getters which have
563 * fewer than AT_MATCH_MINVOTES votes. This may decrease the
564 * number of valid matched pairs below 'nbright', so we
565 * re-set nbright if necessary.
566 */
567 for (i = 0; i < nbright; i++) {
568 if (winner_votes[i] < AT_MATCH_MINVOTES) {
569 #ifdef DEBUG
570 printf("disqualifying all winners after number %d, nbright now %d\n",
571 i, i);
572 #endif
573 nbright = i;
574 break;
575 }
576 }
577
578
579 /*
580 * next, we take the "top" matched pairs of coodinates, and
581 * figure out a transformation of the form
582 *
583 * x' = A + Bx + Cx
584 * y' = D + Ex + Fy
585 *
586 * (i.e. a TRANS structure) which converts the coordinates
587 * of objects in chainA to those in chainB.
588 */
589 if (iter_trans(nbright, star_array_A, num_stars_A,
590 star_array_B, num_stars_B,
591 winner_votes, winner_index_A, winner_index_B,
592 RECALC_NO, max_iter, halt_sigma, trans) != SH_SUCCESS) {
593
594 shError("atFindTrans: iter_trans unable to create a valid TRANS");
595 free_star_array(star_array_A);
596 free_star_array(star_array_B);
597 return(SH_GENERIC_ERROR);
598 }
599
600 #ifdef DEBUG
601 printf(" after calculating new TRANS structure, here it is\n");
602 print_trans(trans);
603 #endif
604
605 /*
606 * clean up memory we allocated during the matching process
607 */
608 free_star_array(star_array_A);
609 free_star_array(star_array_B);
610
611 return(SH_SUCCESS);
612 }
613
614
615
616
617 /************************************************************************
618 * <AUTO EXTRACT>
619 *
620 * ROUTINE: atRecalcTrans
621 *
622 * DESCRIPTION:
623 * Given two lists of stars which ALREADY have been matched,
624 * this routine finds a coord transformation which takes coords
625 * of stars in list A to those in list B.
626 *
627 * We can skip all the matching-triangles business, which makes this
628 * _much_ faster than atFindTrans.
629 *
630 * RETURN:
631 * SH_SUCCESS if all goes well
632 * SH_GENERIC_ERROR if an error occurs
633 *
634 * </AUTO>
635 */
636
637 int
atRecalcTrans(int numA,struct s_star * listA,int numB,struct s_star * listB,int max_iter,double halt_sigma,TRANS * trans)638 atRecalcTrans
639 (
640 int numA, /* I: number of stars in list A */
641 struct s_star *listA, /* I: match this set of objects with list B */
642 int numB, /* I: number of stars in list B */
643 struct s_star *listB, /* I: match this set of objects with list A */
644 int max_iter, /* I: go through at most this many iterations */
645 /* in the iter_trans() loop. */
646 double halt_sigma, /* I: halt the fitting procedure if the mean */
647 /* residual becomes this small */
648 TRANS *trans /* O: place into this TRANS structure's fields */
649 /* the coeffs which convert coords of chainA */
650 /* into coords of chainB system. */
651 )
652 {
653 int i, nbright, min;
654 int num_stars_A; /* number of stars in chain A */
655 int num_stars_B; /* number of stars in chain B */
656 int *winner_votes; /* # votes gotten by top pairs of matched stars */
657 int *winner_index_A; /* elem i in this array is index in star array A */
658 /* which matches ... */
659 int *winner_index_B; /* elem i in this array, index in star array B */
660 int start_pairs;
661 s_star *star_array_A = NULL;
662 s_star *star_array_B = NULL;
663
664 num_stars_A = numA;
665 num_stars_B = numB;
666 star_array_A = list_to_array(numA, listA);
667 star_array_B = list_to_array(numB, listB);
668
669 shAssert(star_array_A != NULL);
670 shAssert(star_array_B != NULL);
671
672 switch (trans->order) {
673 case AT_TRANS_LINEAR:
674 start_pairs = AT_MATCH_STARTN_LINEAR;
675 break;
676 case AT_TRANS_QUADRATIC:
677 start_pairs = AT_MATCH_STARTN_QUADRATIC;
678 break;
679 case AT_TRANS_CUBIC:
680 start_pairs = AT_MATCH_STARTN_CUBIC;
681 break;
682 default:
683 shError("atRecalcTrans: invalid trans->order %d ", trans->order);
684 break;
685 }
686
687 /*
688 * here we check to see if each list of stars contains a
689 * required minimum number of stars. If not, we return with
690 * an error message, and SH_GENERIC_ERROR.
691 *
692 * We set 'nbright' to the minimum of the two list lengths
693 */
694 min = (num_stars_A < num_stars_B ? num_stars_A : num_stars_B);
695 if (min < start_pairs) {
696 shError("atRecalcTrans: only %d stars in list(s), require at least %d",
697 min, start_pairs);
698 free_star_array(star_array_A);
699 free_star_array(star_array_B);
700 return(SH_GENERIC_ERROR);
701 }
702 nbright = min;
703
704 /* this is a sanity check on the above checks */
705 shAssert((nbright >= start_pairs) && (nbright <= min));
706
707
708 #ifdef DEBUG
709 printf("here comes star array A\n");
710 print_star_array(star_array_A, num_stars_A);
711 printf("here comes star array B\n");
712 print_star_array(star_array_B, num_stars_B);
713 #endif
714
715
716 /*
717 * We need to create dummy arrays for 'winner_votes', and the
718 * 'winner_index' arrays. We already know that all these stars
719 * are good matches, and so we can just create some arrays
720 * and fill them with identical numbers. They aren't used by
721 * iter_trans(), anyway.
722 */
723 winner_votes = (int *) shMalloc(nbright*sizeof(int));
724 winner_index_A = (int *) shMalloc(nbright*sizeof(int));
725 winner_index_B = (int *) shMalloc(nbright*sizeof(int));
726 for (i = 0; i < nbright; i++) {
727 winner_votes[i] = 100;
728 winner_index_A[i] = i;
729 winner_index_B[i] = i;
730 }
731
732 /*
733 * next, we take ALL the matched pairs of coodinates, and
734 * figure out a transformation of the form
735 *
736 * x' = A + Bx + Cx
737 * y' = D + Ex + Fy
738 *
739 * (i.e. a TRANS structure) which converts the coordinates
740 * of objects in list A to those in list B
741 */
742 if (iter_trans(nbright, star_array_A, num_stars_A,
743 star_array_B, num_stars_B,
744 winner_votes, winner_index_A, winner_index_B,
745 RECALC_YES, max_iter, halt_sigma, trans) != SH_SUCCESS) {
746
747 shError("atRecalcTrans: iter_trans unable to create a valid TRANS");
748 free_star_array(star_array_A);
749 free_star_array(star_array_B);
750 return(SH_GENERIC_ERROR);
751 }
752
753 #ifdef DEBUG
754 printf(" after calculating new TRANS structure, here it is\n");
755 print_trans(trans);
756 #endif
757
758 /*
759 * clean up memory we allocated during the matching process
760 */
761 shFree(winner_votes);
762 shFree(winner_index_A);
763 shFree(winner_index_B);
764 free_star_array(star_array_A);
765 free_star_array(star_array_B);
766
767 return(SH_SUCCESS);
768 }
769
770
771
772 /************************************************************************
773 * <AUTO EXTRACT>
774 *
775 * ROUTINE: atApplyTrans
776 *
777 * DESCRIPTION:
778 * Given a list of s_star structures, apply the given TRANS to each item in
779 * the list, modifying the "x" and "y" values.
780 *
781 * The TRANS structure has 6 coefficients, which are used as follows:
782 *
783 * x' = A + Bx + Cx
784 * y' = D + Ex + Fy
785 *
786 * Actually, this is a top-level "driver" routine that calls smaller
787 * functions to perform actual tasks. It mostly creates the proper
788 * inputs and outputs for the smaller routines.
789 *
790 *
791 * RETURN:
792 * SH_SUCCESS if all goes well
793 * SH_GENERIC_ERROR if an error occurs
794 *
795 * </AUTO>
796 */
797
798 int
atApplyTrans(int num,s_star * star_list,TRANS * trans)799 atApplyTrans
800 (
801 int num, /* I: number of stars in the linked list */
802 s_star *star_list, /* I/O: modify x,y coords of objects in this list */
803 TRANS *trans /* I: use this TRANS to transform the coords of */
804 /* items in the list */
805 )
806 {
807 int i;
808 struct s_star *ptr, *star_array;
809
810 shAssert(star_list != NULL);
811
812 /*
813 * convert the linked list to an array
814 */
815 star_array = list_to_array(num, star_list);
816
817 #ifdef DEBUG
818 printf("before applying TRANS \n");
819 print_star_array(star_array, num);
820 #endif
821
822 /*
823 * next, apply the transformation to each element of the array
824 */
825 apply_trans(star_array, num, trans);
826
827 #ifdef DEBUG
828 printf("after applying TRANS \n");
829 print_star_array(star_array, num);
830 #endif
831
832 /*
833 * transfer the coord values from the array back into the list
834 */
835 for (ptr = star_list, i = 0; i < num; i++, ptr = ptr->next) {
836 shAssert(ptr != NULL);
837 ptr->x = star_array[i].x;
838 ptr->y = star_array[i].y;
839 }
840
841 /*
842 * delete the array
843 */
844 free_star_array(star_array);
845
846 /*
847 * all done!
848 */
849
850 return(SH_SUCCESS);
851 }
852
853
854 /************************************************************************
855 * <AUTO EXTRACT>
856 *
857 * ROUTINE: atMatchLists
858 *
859 * DESCRIPTION:
860 * Given 2 lists of s_star structures,
861 * which have ALREADY been transformed so that the "x"
862 * and "y" coordinates of each list are close to each other
863 * (i.e. matching items from each list have very similar "x" and "y")
864 * this routine attempts to find all instances of matching items
865 * from the 2 lists.
866 *
867 * We consider a "match" to be the closest coincidence of centers
868 * which are within "radius" pixels of each other.
869 *
870 * Use a slow, but sure, algorithm.
871 *
872 * We will match objects from A --> B. It is possible to have several
873 * As that match to the same B:
874 *
875 * A1 -> B5 and A2 -> B5
876 *
877 * This function finds such multiple-match items and deletes all but
878 * the closest of the matches.
879 *
880 * place the elems of A that are matches into output list J
881 * B that are matches into output list K
882 * A that are not matches into output list L
883 * B that are not matches into output list M
884 *
885 * Place a count of the number of matching pairs into the final
886 * argument, 'num_matches'.
887 *
888 *
889 * RETURN:
890 * SH_SUCCESS if all goes well
891 * SH_GENERIC_ERROR if an error occurs
892 *
893 * </AUTO>
894 */
895
896 int
atMatchLists(int numA,s_star * listA,int numB,s_star * listB,double radius,char * basename,int * num_matches)897 atMatchLists
898 (
899 int numA, /* I: number of stars in list A */
900 s_star *listA, /* I: first input list of items to be matched */
901 int numB, /* I: number of stars in list B */
902 s_star *listB, /* I: second list of items to be matched */
903 double radius, /* I: maximum radius for items to be a match */
904 char *basename, /* I: base of filenames used to store the */
905 /* output; extension indicates which */
906 /* .mtA items from A that matched */
907 /* .mtB items from B that matched */
908 /* .unA items from A that didn't match */
909 /* .unB items from A that didn't match */
910 int *num_matches /* O: number of matching pairs we find */
911 )
912 {
913 s_star *star_array_A;
914 int num_stars_A;
915 s_star *star_array_B;
916 int num_stars_B;
917 s_star *star_array_J, *star_array_K, *star_array_L, *star_array_M;
918 int num_stars_J, num_stars_K, num_stars_L, num_stars_M;
919 char filename[100];
920
921 shAssert(listA != NULL);
922 shAssert(listB != NULL);
923
924 /* convert from linked lists to arrays */
925 num_stars_A = numA;
926 num_stars_B = numB;
927 star_array_A = list_to_array(numA, listA);
928 star_array_B = list_to_array(numB, listB);
929 shAssert(star_array_A != NULL);
930 shAssert(star_array_B != NULL);
931
932 /* reset the 'id' fields in the arrays to match those in the lists */
933 reset_array_ids(listA, numA, star_array_A);
934 reset_array_ids(listB, numB, star_array_B);
935
936
937 /* do the matching process */
938 if (match_arrays_slow(star_array_A, num_stars_A,
939 star_array_B, num_stars_B,
940 radius,
941 &star_array_J, &num_stars_J,
942 &star_array_K, &num_stars_K,
943 &star_array_L, &num_stars_L,
944 &star_array_M, &num_stars_M) != SH_SUCCESS) {
945 shError("atMatchLists: match_arrays_slow fails");
946 return(SH_GENERIC_ERROR);
947 }
948
949 /*
950 * Set the 'num_matches' value to the number of matching pairs
951 * (we could just as well use num_stars_K)
952 */
953 *num_matches = num_stars_J;
954
955 /*
956 * now write the output into ASCII text files, each of which starts
957 * with 'basename', but has a different extension.
958 *
959 * basename.mtA stars from list A that did match array J
960 * basename.mtB stars from list A that did match array K
961 * basename.unA stars from list A that did NOT match array L
962 * basename.unB stars from list A that did NOT match array M
963 */
964 sprintf(filename, "%s.mtA", basename);
965 write_array(num_stars_J, star_array_J, filename);
966 sprintf(filename, "%s.mtB", basename);
967 write_array(num_stars_K, star_array_K, filename);
968 sprintf(filename, "%s.unA", basename);
969 write_array(num_stars_L, star_array_L, filename);
970 sprintf(filename, "%s.unB", basename);
971 write_array(num_stars_M, star_array_M, filename);
972
973 /*
974 * all done!
975 */
976 free_star_array(star_array_J);
977 free_star_array(star_array_K);
978 free_star_array(star_array_L);
979 free_star_array(star_array_M);
980
981 return(SH_SUCCESS);
982 }
983
984
985
986 /************************************************************************
987 * <AUTO EXTRACT>
988 *
989 * ROUTINE: atBuildSmallFile
990 *
991 * DESCRIPTION:
992 * This function is basically the first half of "atFindTrans". It takes
993 * a single list of s_star structures, performs some checks, and calculates
994 * the array of triangles for the bright stars in the list.
995 *
996 * At that point, it creates a new file with name "outfile", and writes
997 * into that file the subset of bright stars and their triangles.
998 * We keep that small file for future use.
999 *
1000 * RETURN:
1001 * SH_SUCCESS if all goes well
1002 * SH_GENERIC_ERROR if an error occurs
1003 *
1004 * </AUTO>
1005 */
1006
1007 int
atBuildSmallFile(double ra,double dec,int numA,struct s_star * listA,int nobj,char * outfile)1008 atBuildSmallFile
1009 (
1010 double ra, /* I: Right Ascension of field center, in degrees */
1011 double dec, /* I: Declination of field center, in degrees */
1012 int numA, /* I: number of stars in list A */
1013 struct s_star *listA, /* I: create an array of triangles for these stars */
1014 int nobj, /* I: max number of bright stars to use in creating */
1015 /* triangles for matching from each list */
1016 char *outfile /* I: create a file with this name, and place lists */
1017 /* of stars and triangles into it */
1018 )
1019 {
1020 int nbright, min;
1021 int num_stars_A; /* number of stars in chain A */
1022 int num_triangles_A; /* number of triangles formed from chain A */
1023 int start_pairs;
1024 s_star *star_array_A = NULL;
1025 s_triangle *triangle_array_A = NULL;
1026
1027 num_stars_A = numA;
1028 star_array_A = list_to_array(numA, listA);
1029
1030 shAssert(star_array_A != NULL);
1031
1032 start_pairs = AT_MATCH_STARTN_LINEAR;
1033
1034 /*
1035 * here we check to see if each list of stars contains a
1036 * required minimum number of stars. If not, we return with
1037 * an error message, and SH_GENERIC_ERROR.
1038 *
1039 * In addition, we check to see that each list has at least 'nobj'
1040 * items. If not, we set 'nbright' to the minimum of the two
1041 * list lengths, and print a warning message so the user knows
1042 * that we're using fewer stars than he asked.
1043 *
1044 * On the other hand, if the user specifies a value of "nobj" which
1045 * is too SMALL, then we ignore it and use the smallest valid
1046 * value (which is start_pairs).
1047 */
1048 min = num_stars_A;
1049 if (min < start_pairs) {
1050 shError("atBuildSmallFile: only %d stars in list, require at least %d",
1051 min, start_pairs);
1052 free_star_array(star_array_A);
1053 return(SH_GENERIC_ERROR);
1054 }
1055 if (nobj > min) {
1056 shDebug(AT_MATCH_ERRLEVEL,
1057 "atBuildSmallFile: using only %d stars, fewer than requested %d",
1058 min, nobj);
1059 nbright = min;
1060 }
1061 else {
1062 nbright = nobj;
1063 }
1064 if (nbright < start_pairs) {
1065 shDebug(AT_MATCH_ERRLEVEL,
1066 "atBuildSmallFile: must use %d stars, more than requested %d",
1067 start_pairs, nobj);
1068 nbright = start_pairs;
1069 }
1070
1071 /* this is a sanity check on the above checks */
1072 shAssert((nbright >= start_pairs) && (nbright <= min));
1073
1074
1075 #ifdef DEBUG
1076 printf("here comes star array A\n");
1077 print_star_array(star_array_A, num_stars_A);
1078 #endif
1079
1080 /*
1081 * we now convert the list of stars into a list of triangles,
1082 * using only a subset of the "nbright" brightest items in the list.
1083 */
1084 triangle_array_A = stars_to_triangles(star_array_A, num_stars_A,
1085 nbright, &num_triangles_A);
1086 shAssert(triangle_array_A != NULL);
1087
1088 /*
1089 * Now we prune the triangle array to eliminate those with
1090 * ratios (b/a) > AT_MATCH_RATIO,
1091 * since Valdes et al. say that this speeds things up and eliminates
1092 * lots of closely-packed triangles.
1093 */
1094 prune_triangle_array(triangle_array_A, &num_triangles_A);
1095 #ifdef DEBUG2
1096 printf("after pruning, here comes triangle array A\n");
1097 print_triangle_array(triangle_array_A, num_triangles_A,
1098 star_array_A, num_stars_A);
1099 #endif
1100
1101 if (write_small_arrays(ra, dec, num_stars_A, star_array_A, nbright,
1102 num_triangles_A, triangle_array_A,
1103 outfile) != SH_SUCCESS) {
1104 free_star_array(star_array_A);
1105 shFree(triangle_array_A);
1106 shError("atBuildSmallFile: write_small_arrays returns with error");
1107 return(SH_GENERIC_ERROR);
1108 }
1109
1110 /* clean up memory */
1111 free_star_array(star_array_A);
1112 shFree(triangle_array_A);
1113
1114 return(SH_SUCCESS);
1115 }
1116
1117
1118 /************************************************************************
1119 * <AUTO EXTRACT>
1120 *
1121 * ROUTINE: atSmallTrans
1122 *
1123 * DESCRIPTION:
1124 * This function is basically the second half of the "atFindTrans"
1125 * function. We pass it a list of detected stars, and a pre-made
1126 * array of catalog stars and triangles.
1127 *
1128 * The first time we call this routine, we convert the _list_ of
1129 * detected stars into an array, and create an array of triangles
1130 * for them. On all subsequent calls, we re-use these arrays.
1131 *
1132 *
1133 * RETURN:
1134 * SH_SUCCESS if all goes well
1135 * SH_GENERIC_ERROR if an error occurs
1136 *
1137 * </AUTO>
1138 */
1139
1140 int
atSmallTrans(int numA,struct s_star * listA,int numB,struct s_star * star_array_B,int num_triangles_B,struct s_triangle * triangle_array_B,double radius,int nobj,double min_scale,double max_scale,double rotation_deg,double tolerance_deg,int max_iter,double halt_sigma,TRANS * trans,int * ntop,int ** top_votes)1141 atSmallTrans
1142 (
1143 int numA, /* I: number of stars in list A */
1144 struct s_star *listA, /* I: match this set of objects with list B */
1145 int numB, /* I: number of stars in pre-made array B */
1146 struct s_star *star_array_B,
1147 /* I: match this array of stars with list A */
1148 int num_triangles_B, /* I: number of pre-made triangles from set B */
1149 struct s_triangle *triangle_array_B,
1150 /* I: pre-made array of triangles */
1151 double radius, /* I: max radius in triangle-space allowed for */
1152 /* a pair of triangles to match */
1153 int nobj, /* I: max num of bright stars to use in creating */
1154 /* triangles for matching from each list */
1155 double min_scale, /* I: minimum permitted relative scale factor */
1156 /* if -1, any scale factor is allowed */
1157 double max_scale, /* I: maximum permitted relative scale factor */
1158 /* if -1, any scale factor is allowed */
1159 double rotation_deg, /* I: desired relative angle of coord systems (deg) */
1160 /* if AT_MATCH_NOANGLE, any orientation is allowed */
1161 double tolerance_deg, /* I: allowed range of orientation angles (deg) */
1162 /* if AT_MATCH_NOANGLE, any orientation is allowed */
1163 int max_iter, /* I: go through at most this many iterations */
1164 /* in the iter_trans() loop. */
1165 double halt_sigma, /* I: halt the fitting procedure if the mean */
1166 /* residual becomes this small */
1167 TRANS *trans, /* O: place into this TRANS structure's fields */
1168 /* the coeffs which convert coords of "A" */
1169 /* into coords of "B" system. */
1170 int *ntop, /* O: number of top "vote getters" in the */
1171 /* top_votes[] array */
1172 int **top_votes /* O: array of votes gotten by the ntop stars */
1173 /* will help evaluate the quality of matches */
1174 )
1175 {
1176 int i, nbright, min, ret;
1177 int num_stars_A; /* number of stars in set A */
1178 int num_stars_B; /* number of stars in set B */
1179 static int num_triangles_A; /* number of triangles formed from set A */
1180 int **vote_matrix;
1181 int *winner_votes; /* # votes gotten by top pairs of matched stars */
1182 int *winner_index_A; /* elem i in this array is index in star array A */
1183 /* which matches ... */
1184 int *winner_index_B; /* elem i in this array, index in star array B */
1185 int start_pairs;
1186 static s_star *star_array_A = NULL;
1187 static s_triangle *triangle_array_A = NULL;
1188 static int first = 1;
1189 int first_flag = 0;
1190
1191 /*
1192 * the first time we call this routine, create an array of stars from
1193 * the items in listA
1194 */
1195 if (first == 1) {
1196 first = 0;
1197 first_flag = 1;
1198 star_array_A = list_to_array(numA, listA);
1199 }
1200
1201 num_stars_A = numA;
1202 num_stars_B = numB;
1203
1204 shAssert(star_array_A != NULL);
1205 shAssert(star_array_B != NULL);
1206
1207 switch (trans->order) {
1208 case AT_TRANS_LINEAR:
1209 start_pairs = AT_MATCH_STARTN_LINEAR;
1210 break;
1211 case AT_TRANS_QUADRATIC:
1212 start_pairs = AT_MATCH_STARTN_QUADRATIC;
1213 break;
1214 case AT_TRANS_CUBIC:
1215 start_pairs = AT_MATCH_STARTN_CUBIC;
1216 break;
1217 default:
1218 shError("atFindTrans: invalid trans->order %d ", trans->order);
1219 break;
1220 }
1221
1222
1223 /*
1224 * here we check to see if each list of stars contains a
1225 * required minimum number of stars. If not, we return with
1226 * an error message, and SH_GENERIC_ERROR.
1227 *
1228 * In addition, we check to see that each list has at least 'nobj'
1229 * items. If not, we set 'nbright' to the minimum of the two
1230 * list lengths, and print a warning message so the user knows
1231 * that we're using fewer stars than he asked.
1232 *
1233 * On the other hand, if the user specifies a value of "nobj" which
1234 * is too SMALL, then we ignore it and use the smallest valid
1235 * value (which is start_pairs).
1236 */
1237 min = (num_stars_A < num_stars_B ? num_stars_A : num_stars_B);
1238 if (min < start_pairs) {
1239 shError("atSmallTrans: only %d stars in list(s), require at least %d",
1240 min, start_pairs);
1241 return(SH_GENERIC_ERROR);
1242 }
1243 if (nobj > min) {
1244 shDebug(AT_MATCH_ERRLEVEL,
1245 "atSmallTrans: using only %d stars, fewer than requested %d",
1246 min, nobj);
1247 nbright = min;
1248 }
1249 else {
1250 nbright = nobj;
1251 }
1252 if (nbright < start_pairs) {
1253 shDebug(AT_MATCH_ERRLEVEL,
1254 "atSmallTrans: must use %d stars, more than requested %d",
1255 start_pairs, nobj);
1256 nbright = start_pairs;
1257 }
1258
1259 /* this is a sanity check on the above checks */
1260 shAssert((nbright >= start_pairs) && (nbright <= min));
1261
1262
1263 #ifdef DEBUG
1264 printf("here comes star array A\n");
1265 print_star_array(star_array_A, num_stars_A);
1266 printf("here comes star array B\n");
1267 print_star_array(star_array_B, num_stars_B);
1268 fflush(stdout);
1269 #endif
1270
1271 /*
1272 * If this is the first time we've entered this function,
1273 * we now convert list A of stars into a list of triangles,
1274 * using only a subset of the "nbright" brightest items in the list.
1275 */
1276 if (first_flag == 1) {
1277 triangle_array_A = stars_to_triangles(star_array_A, num_stars_A,
1278 nbright, &num_triangles_A);
1279 }
1280 shAssert(triangle_array_A != NULL);
1281 shAssert(triangle_array_B != NULL);
1282
1283
1284 /*
1285 * If this is the first time through the function,
1286 * we now prune the "A" triangle arrays to eliminate those with
1287 * ratios (b/a) > AT_MATCH_RATIO,
1288 * since Valdes et al. say that this speeds things up and eliminates
1289 * lots of closely-packed triangles.
1290 *
1291 * Triangle array "B" should have been pruned when it was created,
1292 * so we don't have to do it again.
1293 */
1294 if (first_flag == 1) {
1295 prune_triangle_array(triangle_array_A, &num_triangles_A);
1296 }
1297 #ifdef DEBUG2
1298 printf("after pruning, here comes triangle array A\n");
1299 print_triangle_array(triangle_array_A, num_triangles_A,
1300 star_array_A, num_stars_A);
1301 printf("after pruning, here comes triangle array B\n");
1302 print_triangle_array(triangle_array_B, num_triangles_B,
1303 star_array_B, num_stars_B);
1304 fflush(stdout);
1305 #endif
1306
1307
1308 /*
1309 * Next, we want to try to match triangles in the two arrays.
1310 * What we do is to create a "vote matrix", which is a 2-D array
1311 * with "nbright"-by-"nbright" cells. The cell with
1312 * coords [i][j] holds the number of matched triangles in which
1313 *
1314 * item [i] in star_array_A matches item [j] in star_array_B
1315 *
1316 * We'll use this "vote_matrix" to figure out a first guess
1317 * at the transformation between coord systems.
1318 *
1319 * Note that if there are fewer than "nbright" stars
1320 * in either list, we'll still make the vote_matrix
1321 * contain "nbright"-by-"nbright" cells ...
1322 * there will just be a lot of cells filled with zero.
1323 */
1324 vote_matrix = make_vote_matrix(star_array_A, num_stars_A,
1325 star_array_B, num_stars_B,
1326 triangle_array_A, num_triangles_A,
1327 triangle_array_B, num_triangles_B,
1328 nbright, radius, min_scale, max_scale,
1329 rotation_deg, tolerance_deg);
1330
1331
1332
1333 /*
1334 * having made the vote_matrix, we next need to pick the
1335 * top 'nbright' vote-getters. We call 'top_vote_getters'
1336 * and are given, in its output arguments, pointers to three
1337 * arrays, each of which has 'nbright' elements pertaining
1338 * to a matched pair of STARS:
1339 *
1340 * winner_votes[] number of votes of winners, in descending order
1341 * winner_index_A[] index of star in star_array_A
1342 * winner_index_B[] index of star in star_array_B
1343 *
1344 * Thus, the pair of stars which matched in the largest number
1345 * of triangles will be
1346 *
1347 * star_array_A[winner_index_A[0]] from array A
1348 * star_array_B[winner_index_A[0]] from array B
1349 *
1350 * and the pair of stars which matched in the second-largest number
1351 * of triangles will be
1352 *
1353 * star_array_A[winner_index_A[1]] from array A
1354 * star_array_B[winner_index_A[1]] from array B
1355 *
1356 * and so on.
1357 */
1358 top_vote_getters(vote_matrix, nbright,
1359 &winner_votes, &winner_index_A, &winner_index_B);
1360
1361
1362 /*
1363 * here, we disqualify any of the top vote-getters which have
1364 * fewer than AT_MATCH_MINVOTES votes. This may decrease the
1365 * number of valid matched pairs below 'nbright', so we
1366 * re-set nbright if necessary.
1367 */
1368 for (i = 0; i < nbright; i++) {
1369 if (winner_votes[i] < AT_MATCH_MINVOTES) {
1370 #ifdef DEBUG
1371 printf("disqualifying all winners after number %d, nbright now %d\n",
1372 i, i);
1373 #endif
1374 nbright = i;
1375 break;
1376 }
1377 }
1378
1379 /*
1380 * next, we take the "top" matched pairs of coodinates, and
1381 * figure out a transformation of the form
1382 *
1383 * x' = A + Bx + Cx
1384 * y' = D + Ex + Fy
1385 *
1386 * (i.e. a TRANS structure) which converts the coordinates
1387 * of objects in chainA to those in chainB.
1388 */
1389 ret = iter_trans(nbright, star_array_A, num_stars_A,
1390 star_array_B, num_stars_B,
1391 winner_votes, winner_index_A, winner_index_B,
1392 RECALC_NO, max_iter, halt_sigma, trans);
1393 if (ret != SH_SUCCESS) {
1394 shDebug(AT_MATCH_ERRLEVEL,
1395 "atSmallTrans: iter_trans unable to create a valid TRANS");
1396 return(SH_GENERIC_ERROR);
1397 }
1398
1399 #ifdef DEBUG
1400 printf(" after calculating new TRANS structure, here it is\n");
1401 print_trans(trans);
1402 #endif
1403
1404 /*
1405 * set the output args "ntop" and "winner_votes"
1406 */
1407 *ntop = nbright;
1408 *top_votes = winner_votes;
1409
1410 return(SH_SUCCESS);
1411 }
1412
1413
1414
1415
1416 /* end of PUBLIC information */
1417 /*-----------------------------------------------------------------------*/
1418 /* start of PRIVATE information */
1419
1420 /*
1421 * the functions listed from here on are intended to be used only
1422 * "internally", called by the PUBLIC functions above. Users
1423 * should be discouraged from accessing them directly.
1424 */
1425
1426
1427 /************************************************************************
1428 *
1429 *
1430 * ROUTINE: set_star
1431 *
1432 * DESCRIPTION:
1433 * Given a pointer to an EXISTING s_star, initialize its values
1434 * and set x, y, and mag to the given values.
1435 *
1436 * RETURN:
1437 * SH_SUCCESS if all goes well
1438 * SH_GENERIC_ERROR if not
1439 *
1440 * </AUTO>
1441 */
1442
1443 static int
set_star(s_star * star,double x,double y,double mag)1444 set_star
1445 (
1446 s_star *star, /* I: pointer to existing s_star structure */
1447 double x, /* I: star's "X" coordinate */
1448 double y, /* I: star's "Y" coordinate */
1449 double mag /* I: star's "mag" coordinate */
1450 )
1451 {
1452 static int id_number = 0;
1453
1454 if (star == NULL) {
1455 shError("set_star: given a NULL star");
1456 return(SH_GENERIC_ERROR);
1457 }
1458 star->id = id_number++;
1459 star->index = -1;
1460 star->x = x;
1461 star->y = y;
1462 star->mag = mag;
1463 star->match_id = -1;
1464 star->next = (s_star *) NULL;
1465
1466 return(SH_SUCCESS);
1467 }
1468
1469
1470
1471 /************************************************************************
1472 *
1473 *
1474 * ROUTINE: copy_star
1475 *
1476 * DESCRIPTION:
1477 * Copy the contents of the "s_star" to which "from_ptr" points
1478 * to the "s_star" to which "to_ptr" points.
1479 *
1480 * RETURN:
1481 * nothing
1482 *
1483 * </AUTO>
1484 */
1485
1486 static void
copy_star(s_star * from_ptr,s_star * to_ptr)1487 copy_star
1488 (
1489 s_star *from_ptr, /* I: copy contents of _this_ star ... */
1490 s_star *to_ptr /* O: into _this_ star */
1491 )
1492 {
1493 shAssert(from_ptr != NULL);
1494 shAssert(to_ptr != NULL);
1495
1496 to_ptr->id = from_ptr->id;
1497 to_ptr->index = from_ptr->index;
1498 to_ptr->x = from_ptr->x;
1499 to_ptr->y = from_ptr->y;
1500 to_ptr->mag = from_ptr->mag;
1501 to_ptr->match_id = from_ptr->match_id;
1502 to_ptr->next = from_ptr->next;
1503
1504 }
1505
1506
1507
1508
1509 /************************************************************************
1510 *
1511 *
1512 * ROUTINE: copy_star_array
1513 *
1514 * DESCRIPTION:
1515 * Given to arrays of "s_star" structures, EACH OF WHICH MUST
1516 * ALREADY HAVE BEEN ALLOCATED and have "num" elements,
1517 * copy the contents of the items in "from_array"
1518 * to those in "to_array".
1519 *
1520 * RETURN:
1521 * nothing
1522 *
1523 * </AUTO>
1524 */
1525
1526 static void
copy_star_array(s_star * from_array,s_star * to_array,int num_stars)1527 copy_star_array
1528 (
1529 s_star *from_array, /* I: copy contents of _this_ array ... */
1530 s_star *to_array, /* O: into _this_ array */
1531 int num_stars /* I: each aray must have this many elements */
1532 )
1533 {
1534 int i;
1535 s_star *from_ptr, *to_ptr;
1536
1537 shAssert(from_array != NULL);
1538 shAssert(to_array != NULL);
1539
1540 for (i = 0; i < num_stars; i++) {
1541 from_ptr = &(from_array[i]);
1542 to_ptr = &(to_array[i]);
1543 shAssert(from_ptr != NULL);
1544 shAssert(to_ptr != NULL);
1545
1546 to_ptr->id = from_ptr->id;
1547 to_ptr->index = from_ptr->index;
1548 to_ptr->x = from_ptr->x;
1549 to_ptr->y = from_ptr->y;
1550 to_ptr->mag = from_ptr->mag;
1551 to_ptr->match_id = from_ptr->match_id;
1552 to_ptr->next = from_ptr->next;
1553 }
1554
1555 }
1556
1557
1558
1559 /************************************************************************
1560 *
1561 *
1562 * ROUTINE: free_star_array
1563 *
1564 * DESCRIPTION:
1565 * Delete an array of "num" s_star structures.
1566 *
1567 * RETURN:
1568 * nothing
1569 *
1570 * </AUTO>
1571 */
1572
1573 static void
free_star_array(s_star * first)1574 free_star_array
1575 (
1576 s_star *first /* first star in the array to be deleted */
1577 )
1578 {
1579 shFree(first);
1580 }
1581
1582
1583
1584
1585 /************************************************************************
1586 *
1587 *
1588 * ROUTINE: print_star_array
1589 *
1590 * DESCRIPTION:
1591 * Given an array of "num" s_star structures, print out
1592 * a bit of information on each in a single line.
1593 *
1594 * For debugging purposes.
1595 *
1596 * RETURN:
1597 * nothing
1598 *
1599 * </AUTO>
1600 */
1601
1602 #ifdef DEBUG
1603
1604 static void
print_star_array(s_star * array,int num)1605 print_star_array
1606 (
1607 s_star *array, /* I: first star in array */
1608 int num /* I: number of stars in the array to print */
1609 )
1610 {
1611 int i;
1612 s_star *star;
1613
1614 for (i = 0; i < num; i++) {
1615 star = &(array[i]);
1616 shAssert(star != NULL);
1617 printf(" %4d %4d %11.4e %11.4e %6.2f\n", i, star->id, star->x,
1618 star->y, star->mag);
1619 }
1620 }
1621
1622 #endif /* DEBUG */
1623
1624
1625 /************************************************************************
1626 *
1627 *
1628 * ROUTINE: calc_distances
1629 *
1630 * DESCRIPTION:
1631 * Given an array of N='numstars' s_star structures, create a 2-D array
1632 * called "matrix" with NxN elements and fill it by setting
1633 *
1634 * matrix[i][j] = distance between stars i and j
1635 *
1636 * where 'i' and 'j' are the indices of their respective stars in
1637 * the 1-D array.
1638 *
1639 * RETURN:
1640 * double **array pointer to array of pointers to each row of array
1641 * NULL if something goes wrong.
1642 *
1643 * </AUTO>
1644 */
1645
1646 static double **
calc_distances(s_star * star_array,int numstars)1647 calc_distances
1648 (
1649 s_star *star_array, /* I: array of s_stars */
1650 int numstars /* I: with this many elements */
1651 )
1652 {
1653 int i, j;
1654 double **matrix;
1655 double dx, dy, dist;
1656
1657 if (numstars == 0) {
1658 shError("calc_distances: given an array of zero stars");
1659 return(NULL);
1660 }
1661
1662 /* allocate the array, row-by-row */
1663 matrix = (double **) shMalloc(numstars*sizeof(double *));
1664 for (i = 0; i < numstars; i++) {
1665 matrix[i] = (double *) shMalloc(numstars*sizeof(double));
1666 }
1667
1668 /* fill up the array */
1669 for (i = 0; i < numstars - 1; i++) {
1670 for (j = i + 1; j < numstars; j++) {
1671 dx = star_array[i].x - star_array[j].x;
1672 dy = star_array[i].y - star_array[j].y;
1673 dist = sqrt(dx*dx + dy*dy);
1674 matrix[i][j] = (double) dist;
1675 matrix[j][i] = (double) dist;
1676 }
1677 }
1678 /* for safety's sake, let's fill the diagonal elements with zeros */
1679 for (i = 0; i < numstars; i++) {
1680 matrix[i][i] = 0.0;
1681 }
1682
1683 /* okay, we're done. return a pointer to the array */
1684 return(matrix);
1685 }
1686
1687
1688 /************************************************************************
1689 *
1690 *
1691 * ROUTINE: free_distances
1692 *
1693 * DESCRIPTION:
1694 * Given a 2-D array of "num"-by-"num" double elements, free up
1695 * each row of the array, then free the array itself.
1696 *
1697 * RETURN:
1698 * nothing
1699 *
1700 * </AUTO>
1701 */
1702
1703 static void
free_distances(double ** array,int num)1704 free_distances
1705 (
1706 double **array, /* I: square array we'll free */
1707 int num /* I: number of elems in each row */
1708 )
1709 {
1710 int i;
1711
1712 for (i = 0; i < num; i++) {
1713 shFree(array[i]);
1714 }
1715 shFree(array);
1716 }
1717
1718
1719 /************************************************************************
1720 *
1721 *
1722 * ROUTINE: print_dist_matrix
1723 *
1724 * DESCRIPTION:
1725 * Given a 2-D array of "num"-by-"num" distances between pairs of
1726 * stars, print out the 2-D array in a neat fashion.
1727 *
1728 * For debugging purposes.
1729 *
1730 * RETURN:
1731 * nothing
1732 *
1733 * </AUTO>
1734 */
1735
1736 #ifdef DEBUG
1737
1738 static void
print_dist_matrix(double ** matrix,int num)1739 print_dist_matrix
1740 (
1741 double **matrix, /* I: pointer to start of 2-D square array */
1742 int num /* I: number of rows and columns in the array */
1743 )
1744 {
1745 int i, j;
1746
1747 for (i = 0; i < num; i++) {
1748 shAssert(matrix[i] != NULL);
1749 for (j = 0; j < num; j++) {
1750 printf("%11.4e ", matrix[i][j]);
1751 }
1752 printf("\n");
1753 }
1754 }
1755
1756 #endif /* DEBUG */
1757
1758
1759
1760 /************************************************************************
1761 *
1762 *
1763 * ROUTINE: set_triangle
1764 *
1765 * DESCRIPTION:
1766 * Set the elements of some given, EXISTING instance of an "s_triangle"
1767 * structure, given (the indices to) three s_star structures for its vertices.
1768 * We check to make sure
1769 * that the three stars are three DIFFERENT stars, asserting
1770 * if not.
1771 *
1772 * The triangle's "a_index" is set to the position of the star opposite
1773 * its side "a" in its star array, and similarly for "b_index" and "c_index".
1774 *
1775 * RETURN:
1776 * nothing
1777 *
1778 * </AUTO>
1779 */
1780
1781 static void
set_triangle(s_triangle * tri,s_star * star_array,int s1,int s2,int s3,double ** darray)1782 set_triangle
1783 (
1784 s_triangle *tri, /* we set fields of this existing structure */
1785 s_star *star_array, /* use stars in this array as vertices */
1786 int s1, /* index in 'star_array' of one vertex */
1787 int s2, /* index in 'star_array' of one vertex */
1788 int s3, /* index in 'star_array' of one vertex */
1789 double **darray /* array of distances between stars */
1790 )
1791 {
1792 static int id_number = 0;
1793 double d12, d23, d13;
1794 double a, b, c;
1795 s_star *star1, *star2, *star3;
1796
1797 shAssert(tri != NULL);
1798 shAssert((s1 != s2) && (s1 != s3) && (s2 != s3));
1799 star1 = &star_array[s1];
1800 star2 = &star_array[s2];
1801 star3 = &star_array[s3];
1802 shAssert((star1 != NULL) && (star2 != NULL) && (star3 != NULL));
1803
1804 tri->id = id_number++;
1805 tri->index = -1;
1806
1807 /*
1808 * figure out which sides is longest and shortest, and assign
1809 *
1810 * "a" to the length of the longest side
1811 * "b" intermediate
1812 * "c" shortest
1813 *
1814 * We use temp variables d12 = distance between stars 1 and 2
1815 * d23 = distance between stars 2 and 3
1816 * d13 = distance between stars 1 and 3
1817 * for convenience.
1818 *
1819 */
1820 d12 = darray[s1][s2];
1821 d23 = darray[s2][s3];
1822 d13 = darray[s1][s3];
1823
1824 /* sanity check */
1825 shAssert(d12 >= 0.0);
1826 shAssert(d23 >= 0.0);
1827 shAssert(d13 >= 0.0);
1828
1829 if ((d12 >= d23) && (d12 >= d13)) {
1830 /* this applies if the longest side connects stars 1 and 2 */
1831 tri->a_index = star3->index;
1832 a = d12;
1833 if (d23 >= d13) {
1834 tri->b_index = star1->index;
1835 b = d23;
1836 tri->c_index = star2->index;
1837 c = d13;
1838 }
1839 else {
1840 tri->b_index = star2->index;
1841 b = d13;
1842 tri->c_index = star1->index;
1843 c = d23;
1844 }
1845 }
1846 else if ((d23 > d12) && (d23 >= d13)) {
1847 /* this applies if the longest side connects stars 2 and 3 */
1848 tri->a_index = star1->index;
1849 a = d23;
1850 if (d12 > d13) {
1851 tri->b_index = star3->index;
1852 b = d12;
1853 tri->c_index = star2->index;
1854 c = d13;
1855 }
1856 else {
1857 tri->b_index = star2->index;
1858 b = d13;
1859 tri->c_index = star3->index;
1860 c = d12;
1861 }
1862 }
1863 else if ((d13 > d12) && (d13 > d23)) {
1864 /* this applies if the longest side connects stars 1 and 3 */
1865 tri->a_index = star2->index;
1866 a = d13;
1867 if (d12 > d23) {
1868 tri->b_index = star3->index;
1869 b = d12;
1870 tri->c_index = star1->index;
1871 c = d23;
1872 }
1873 else {
1874 tri->b_index = star1->index;
1875 b = d23;
1876 tri->c_index = star3->index;
1877 c = d12;
1878 }
1879 }
1880 else {
1881 /* we should never get here! */
1882 shError("set_triangle: impossible situation?!");
1883 shAssert(0);
1884 }
1885
1886 /*
1887 * now that we've figured out the longest, etc., sides, we can
1888 * fill in the rest of the triangle's elements
1889 *
1890 * We need to make a special check, in case a == 0. In that
1891 * case, we'll just set the ratios ba and ca = 1.0, and hope
1892 * that these triangles are ignored.
1893 */
1894 tri->a_length = a;
1895 if (a > 0.0) {
1896 tri->ba = b/a;
1897 tri->ca = c/a;
1898 }
1899 else {
1900 tri->ba = 1.0;
1901 tri->ca = 1.0;
1902 }
1903 tri->side_a_angle = atan2(star_array[tri->a_index].y - star_array[tri->b_index].y,
1904 star_array[tri->a_index].x - star_array[tri->b_index].x);
1905 #ifdef DEBUG2
1906 printf(" triangle %5d has side_a_angle %lf = %lf deg \n",
1907 tri->id, tri->side_a_angle, tri->side_a_angle*(180/3.14159));
1908 #endif
1909
1910 tri->match_id = -1;
1911 tri->next = (s_triangle *) NULL;
1912 }
1913
1914
1915 /************************************************************************
1916 *
1917 *
1918 * ROUTINE: print_triangle_array
1919 *
1920 * DESCRIPTION:
1921 * Given an array of "numtriangle" s_triangle structures,
1922 * and an array of "numstars" s_star structures that make them up,
1923 * print out
1924 * a bit of information on each triangle in a single line.
1925 *
1926 * For debugging purposes.
1927 *
1928 * RETURN:
1929 * nothing
1930 *
1931 * </AUTO>
1932 */
1933
1934 #ifdef DEBUG2
1935
1936 static void
print_triangle_array(s_triangle * t_array,int numtriangles,s_star * star_array,int numstars)1937 print_triangle_array
1938 (
1939 s_triangle *t_array, /* I: first triangle in array */
1940 int numtriangles, /* I: number of triangles in the array to print */
1941 s_star *star_array, /* I: array of stars which appear in triangles */
1942 int numstars /* I: number of stars in star_array */
1943 )
1944 {
1945 int i;
1946 s_triangle *triangle;
1947 s_star *sa, *sb, *sc;
1948
1949 for (i = 0; i < numtriangles; i++) {
1950 triangle = &(t_array[i]);
1951 shAssert(triangle != NULL);
1952
1953 sa = &(star_array[triangle->a_index]);
1954 sb = &(star_array[triangle->b_index]);
1955 sc = &(star_array[triangle->c_index]);
1956
1957 printf("%4d %4d %3d (%5.1f,%5.1f) %3d (%5.1f,%5.1f) %3d (%5.1f, %5.1f) %5.3f %5.3f\n",
1958 i, triangle->id,
1959 triangle->a_index, sa->x, sa->y,
1960 triangle->b_index, sb->x, sb->y,
1961 triangle->c_index, sc->x, sc->y,
1962 triangle->ba, triangle->ca);
1963 }
1964 }
1965
1966 #endif /* DEBUG */
1967
1968
1969
1970
1971 /************************************************************************
1972 *
1973 *
1974 * ROUTINE: stars_to_triangles
1975 *
1976 * DESCRIPTION:
1977 * Convert an array of s_stars to an array of s_triangles.
1978 * We use only the brightest 'nbright' objects in the linked list.
1979 * The steps we need to take are:
1980 *
1981 * 1. sort the array of s_stars by magnitude, and
1982 * set "index" values in the sorted list.
1983 * 2. calculate star-star distances in the sorted list,
1984 * (for the first 'nbright' objects only)
1985 * (creates a 2-D array of distances)
1986 * 3. create a linked list of all possible triangles
1987 * (again using the first 'nbright' objects only)
1988 * 4. clean up -- delete the 2-D array of distances
1989 *
1990 * We place the number of triangles in the final argument, and
1991 * return a pointer to the new array of s_triangle structures.
1992 *
1993 * RETURN:
1994 * s_triangle * pointer to new array of triangles
1995 * (and # of triangles put into output arg)
1996 * NULL if error occurs
1997 *
1998 * </AUTO>
1999 */
2000
2001 static s_triangle *
stars_to_triangles(s_star * star_array,int numstars,int nbright,int * numtriangles)2002 stars_to_triangles
2003 (
2004 s_star *star_array, /* I: array of s_stars */
2005 int numstars, /* I: the total number of stars in the array */
2006 int nbright, /* I: use only the 'nbright' brightest stars */
2007 int *numtriangles /* O: number of triangles we create */
2008 )
2009 {
2010 int numt;
2011 double **dist_matrix;
2012 s_triangle *triangle_array;
2013
2014 /*
2015 * check to see if 'nbright' > 'numstars' ... if so, we re-set
2016 * nbright = numstars
2017 *
2018 * so that we don't have to try to keep track of them separately.
2019 * We'll be able to use 'nbright' safely from then on in this function.
2020 */
2021 if (numstars < nbright) {
2022 nbright = numstars;
2023 }
2024
2025 /*
2026 * sort the stars in the array by their 'mag' field, so that we get
2027 * them in order "brightest-first".
2028 */
2029 sort_star_by_mag(star_array, numstars);
2030
2031 #ifdef DEBUG
2032 printf("stars_to_triangles: here comes star array after sorting\n");
2033 print_star_array(star_array, numstars);
2034 #endif
2035
2036 /*
2037 * calculate the distances between each pair of stars, placing them
2038 * into the newly-created 2D array called "dist_matrix". Note that
2039 * we only need to include the first 'nbright' stars in the
2040 * distance calculations.
2041 */
2042 dist_matrix = calc_distances(star_array, nbright);
2043 shAssert(dist_matrix != NULL);
2044
2045 #ifdef DEBUG
2046 printf("stars_to_triangles: here comes distance matrix\n");
2047 print_dist_matrix(dist_matrix, nbright);
2048 #endif
2049
2050
2051 /*
2052 * create an array of the appropriate number of triangles that
2053 * can be formed from the 'nbright' objects.
2054 */
2055 numt = (nbright*(nbright - 1)*(nbright - 2))/6;
2056 *numtriangles = numt;
2057 triangle_array = (s_triangle *) shMalloc(numt*sizeof(s_triangle));
2058
2059
2060 /*
2061 * now let's fill that array by making all the possible triangles
2062 * out of the first 'nbright' objects in the array of stars.
2063 */
2064 fill_triangle_array(star_array, nbright, dist_matrix,
2065 *numtriangles, triangle_array);
2066
2067 #ifdef DEBUG2
2068 printf("stars_to_triangles: here comes the triangle array\n");
2069 print_triangle_array(triangle_array, *numtriangles, star_array, nbright);
2070 #endif
2071
2072
2073 /*
2074 * we've successfully created the array of triangles, so we can
2075 * now get rid of the "dist_matrix" array. We won't need it
2076 * any more.
2077 */
2078 free_distances(dist_matrix, nbright);
2079
2080 return(triangle_array);
2081 }
2082
2083
2084 /************************************************************************
2085 *
2086 *
2087 * ROUTINE: sort_star_by_mag
2088 *
2089 * DESCRIPTION:
2090 * Given an array of "num" s_star structures, sort it in order
2091 * of increasing magnitude.
2092 *
2093 * After sorting, walk through the array and set each star's
2094 * "index" field equal to the star's position in the array.
2095 * Thus, the first star will have index=0, and the second index=1,
2096 * and so forth.
2097 *
2098 * Calls the "compare_star_by_mag" function, below.
2099 *
2100 * RETURN:
2101 * nothing
2102 *
2103 * </AUTO>
2104 */
2105
2106 static void
sort_star_by_mag(s_star * array,int num)2107 sort_star_by_mag
2108 (
2109 s_star *array, /* I: array of structures to be sorted */
2110 int num /* I: number of stars in the array */
2111 )
2112 {
2113 int i;
2114
2115 qsort((char *) array, num, sizeof(s_star), (PFI) compare_star_by_mag);
2116
2117 /* now set the "index" field for each star */
2118 for (i = 0; i < num; i++) {
2119 array[i].index = i;
2120 }
2121 }
2122
2123 /************************************************************************
2124 *
2125 *
2126 * ROUTINE: compare_star_by_mag
2127 *
2128 * DESCRIPTION:
2129 * Given two s_star structures, compare their "mag" values.
2130 * Used by "sort_star_by_mag".
2131 *
2132 * RETURN:
2133 * 1 if first star has larger "mag"
2134 * 0 if the two have equal "mag"
2135 * -1 if the first has smaller "mag"
2136 *
2137 * </AUTO>
2138 */
2139
2140 static int
compare_star_by_mag(s_star * star1,s_star * star2)2141 compare_star_by_mag
2142 (
2143 s_star *star1, /* I: compare "mag" field of THIS star ... */
2144 s_star *star2 /* I: ... with THIS star */
2145 )
2146 {
2147 shAssert((star1 != NULL) && (star2 != NULL));
2148
2149 if (star1->mag > star2->mag) {
2150 return(1);
2151 }
2152 if (star1->mag < star2->mag) {
2153 return(-1);
2154 }
2155 return(0);
2156 }
2157
2158
2159 /************************************************************************
2160 *
2161 *
2162 * ROUTINE: sort_star_by_x
2163 *
2164 * DESCRIPTION:
2165 * Given an array of "num" s_star structures, sort it in order
2166 * of increasing "x" values.
2167 *
2168 * In this case, we do NOT re-set the "index" field of each
2169 * s_star after sorting!
2170 *
2171 * Calls the "compare_star_by_x" function, below.
2172 *
2173 * RETURN:
2174 * nothing
2175 *
2176 * </AUTO>
2177 */
2178
2179 static void
sort_star_by_x(s_star * array,int num)2180 sort_star_by_x
2181 (
2182 s_star *array, /* I: array of structures to be sorted */
2183 int num /* I: number of stars in the array */
2184 )
2185 {
2186 qsort((char *) array, num, sizeof(s_star), (PFI) compare_star_by_x);
2187 }
2188
2189 /************************************************************************
2190 *
2191 *
2192 * ROUTINE: compare_star_by_x
2193 *
2194 * DESCRIPTION:
2195 * Given two s_star structures, compare their "x" values.
2196 * Used by "sort_star_by_x".
2197 *
2198 * RETURN:
2199 * 1 if first star has larger "x"
2200 * 0 if the two have equal "x"
2201 * -1 if the first has smaller "x"
2202 *
2203 * </AUTO>
2204 */
2205
2206 static int
compare_star_by_x(s_star * star1,s_star * star2)2207 compare_star_by_x
2208 (
2209 s_star *star1, /* I: compare "x" field of THIS star ... */
2210 s_star *star2 /* I: ... with THIS star */
2211 )
2212 {
2213 shAssert((star1 != NULL) && (star2 != NULL));
2214
2215 if (star1->x > star2->x) {
2216 return(1);
2217 }
2218 if (star1->x < star2->x) {
2219 return(-1);
2220 }
2221 return(0);
2222 }
2223
2224 /************************************************************************
2225 *
2226 *
2227 * ROUTINE: sort_star_by_match_id
2228 *
2229 * DESCRIPTION:
2230 * Given an array of "num" s_star structures, sort it in order
2231 * of increasing "match_id" values.
2232 *
2233 * In this case, we do NOT re-set the "index" field of each
2234 * s_star after sorting!
2235 *
2236 * Calls the "compare_star_by_match_id" function, below.
2237 *
2238 * RETURN:
2239 * nothing
2240 *
2241 * </AUTO>
2242 */
2243
2244 static void
sort_star_by_match_id(s_star * array,int num)2245 sort_star_by_match_id
2246 (
2247 s_star *array, /* I: array of structures to be sorted */
2248 int num /* I: number of stars in the array */
2249 )
2250 {
2251 qsort((char *) array, num, sizeof(s_star), (PFI) compare_star_by_match_id);
2252 }
2253
2254 /************************************************************************
2255 *
2256 *
2257 * ROUTINE: compare_star_by_match_id
2258 *
2259 * DESCRIPTION:
2260 * Given two s_star structures, compare their "match_id" values.
2261 * Used by "sort_star_by_match_id".
2262 *
2263 * RETURN:
2264 * 1 if first star has larger "match_id"
2265 * 0 if the two have equal "match_id"
2266 * -1 if the first has smaller "match_id"
2267 *
2268 * </AUTO>
2269 */
2270
2271 static int
compare_star_by_match_id(s_star * star1,s_star * star2)2272 compare_star_by_match_id
2273 (
2274 s_star *star1, /* I: compare "match_id" field of THIS star ... */
2275 s_star *star2 /* I: ... with THIS star */
2276 )
2277 {
2278 shAssert((star1 != NULL) && (star2 != NULL));
2279
2280 if (star1->match_id > star2->match_id) {
2281 return(1);
2282 }
2283 if (star1->match_id < star2->match_id) {
2284 return(-1);
2285 }
2286 return(0);
2287 }
2288
2289
2290
2291 /************************************************************************
2292 *
2293 *
2294 * ROUTINE: fill_triangle_array
2295 *
2296 * DESCRIPTION:
2297 * Given an array of stars, and a matrix of distances between them,
2298 * form all the triangles possible; place the properties of these
2299 * triangles into the "t_array" array, which must already have
2300 * been allocated and contain "numtriangles" elements.
2301 *
2302 * RETURN:
2303 * SH_SUCCESS if all goes well
2304 * SH_GENERIC_ERROR if error occurs
2305 *
2306 * </AUTO>
2307 */
2308
2309 static int
fill_triangle_array(s_star * star_array,int numstars,double ** dist_matrix,int numtriangles,s_triangle * t_array)2310 fill_triangle_array
2311 (
2312 s_star *star_array, /* I: array of stars we use to form triangles */
2313 int numstars, /* I: use this many stars from the array */
2314 double **dist_matrix, /* I: numstars-by-numstars matrix of distances */
2315 /* between stars in the star_array */
2316 int numtriangles, /* I: number of triangles in the t_array */
2317 s_triangle *t_array /* O: we'll fill properties of triangles in */
2318 /* this array, which must already exist */
2319 )
2320 {
2321 int i, j, k, n;
2322 s_triangle *triangle;
2323
2324 shAssert((star_array != NULL) && (dist_matrix != NULL) && (t_array != NULL));
2325
2326
2327 n = 0;
2328 for (i = 0; i < numstars - 2; i++) {
2329 for (j = i + 1; j < numstars - 1; j++) {
2330 for (k = j + 1; k < numstars; k++) {
2331
2332 triangle = &(t_array[n]);
2333 set_triangle(triangle, star_array, i, j, k, dist_matrix);
2334
2335 n++;
2336 }
2337 }
2338 }
2339 shAssert(n == numtriangles);
2340
2341 return(SH_SUCCESS);
2342 }
2343
2344
2345 /************************************************************************
2346 *
2347 *
2348 * ROUTINE: sort_triangle_array
2349 *
2350 * DESCRIPTION:
2351 * Given an array of "num" s_triangle structures, sort it in order
2352 * of increasing "ba" value (where "ba" is the ratio of lengths
2353 * of side b to side a).
2354 *
2355 * Calls the "compare_triangle" function, below.
2356 *
2357 * RETURN:
2358 * nothing
2359 *
2360 * </AUTO>
2361 */
2362
2363 static void
sort_triangle_array(s_triangle * array,int num)2364 sort_triangle_array
2365 (
2366 s_triangle *array, /* I: array of structures to be sorted */
2367 int num /* I: number of triangles in the array */
2368 )
2369 {
2370 qsort((char *) array, num, sizeof(s_triangle), (PFI) compare_triangle);
2371 }
2372
2373
2374 /************************************************************************
2375 *
2376 *
2377 * ROUTINE: compare_triangle
2378 *
2379 * DESCRIPTION:
2380 * Given two s_triangle structures, compare their "ba" values.
2381 * Used by "sort_triangle_array".
2382 *
2383 * RETURN:
2384 * 1 if first star has larger "ba"
2385 * 0 if the two have equal "ba"
2386 * -1 if the first has smaller "ba"
2387 *
2388 * </AUTO>
2389 */
2390
2391 static int
compare_triangle(s_triangle * triangle1,s_triangle * triangle2)2392 compare_triangle
2393 (
2394 s_triangle *triangle1, /* I: compare "ba" field of THIS triangle ... */
2395 s_triangle *triangle2 /* I: ... with THIS triangle */
2396 )
2397 {
2398 shAssert((triangle1 != NULL) && (triangle2 != NULL));
2399
2400 if (triangle1->ba > triangle2->ba) {
2401 return(1);
2402 }
2403 if (triangle1->ba < triangle2->ba) {
2404 return(-1);
2405 }
2406 return(0);
2407 }
2408
2409 /************************************************************************
2410 *
2411 *
2412 * ROUTINE: find_ba_triangle
2413 *
2414 * DESCRIPTION:
2415 * Given an array of "num" s_triangle structures, which have already
2416 * been sorted in order of increasing "ba" value, and given one
2417 * particular "ba" value ba0, return the index of the first triangle
2418 * in the array which has "ba" >= ba0.
2419 *
2420 * We use a binary search, on the "ba" element of each structure.
2421 *
2422 * If there is no such triangle, just return the index of the last
2423 * triangle in the list.
2424 *
2425 * Calls the "compare_triangle" function, above.
2426 *
2427 * RETURN:
2428 * index of closest triangle in array if all goes well
2429 * index of last triangle in array if nothing close
2430 *
2431 * </AUTO>
2432 */
2433
2434 static int
find_ba_triangle(s_triangle * array,int num,double ba0)2435 find_ba_triangle
2436 (
2437 s_triangle *array, /* I: array of structures which been sorted */
2438 int num, /* I: number of triangles in the array */
2439 double ba0 /* I: value of "ba" we seek */
2440 )
2441 {
2442 int top, bottom, mid;
2443
2444 #ifdef DEBUG2
2445 printf("find_ba_triangle: looking for ba = %.2f\n", ba0);
2446 #endif
2447
2448 top = 0;
2449 if ((bottom = num - 1) < 0) {
2450 bottom = 0;
2451 }
2452
2453 while (bottom - top > 2) {
2454 mid = (top + bottom)/2;
2455 #ifdef DEBUG2
2456 printf(" array[%4d] ba=%.2f array[%4d] ba=%.2f array[%4d] ba=%.2f\n",
2457 top, array[top].ba, mid, array[mid].ba,
2458 bottom, array[bottom].ba);
2459 #endif
2460 if (array[mid].ba < ba0) {
2461 top = mid;
2462 }
2463 else {
2464 bottom = mid;
2465 }
2466 }
2467 #ifdef DEBUG2
2468 printf(" array[%4d] ba=%.2f array[%4d] ba=%.2f\n",
2469 top, array[top].ba, bottom, array[bottom].ba);
2470 #endif
2471
2472 /*
2473 * if we get here, then the item we seek is either "top" or "bottom"
2474 * (which may point to the same item in the array).
2475 */
2476 if (array[top].ba < ba0) {
2477 #ifdef DEBUG2
2478 printf(" returning array[%4d] ba=%.2f \n", bottom, array[bottom].ba);
2479 #endif
2480 return(bottom);
2481 }
2482 else {
2483 #ifdef DEBUG2
2484 printf(" returning array[%4d] ba=%.2f \n", top, array[top].ba);
2485 #endif
2486 return(top);
2487 }
2488 }
2489
2490
2491
2492
2493 /************************************************************************
2494 *
2495 *
2496 * ROUTINE: prune_triangle_array
2497 *
2498 * DESCRIPTION:
2499 * Given an array of triangles, sort them in increasing order
2500 * of the side ratio (b/a), and then "ignore" all triangles
2501 * with (b/a) > AT_MATCH_RATIO.
2502 *
2503 * We re-set the arg "numtriangles" as needed, but leave the
2504 * space in the array allocated (since the array was allocated
2505 * as a single block).
2506 *
2507 * RETURN:
2508 * nothing
2509 *
2510 * </AUTO>
2511 */
2512
2513 static void
prune_triangle_array(s_triangle * t_array,int * numtriangles)2514 prune_triangle_array
2515 (
2516 s_triangle *t_array, /* I/O: array of triangles to sort and prune */
2517 int *numtriangles /* I/O: number of triangles in the t_array */
2518 )
2519 {
2520 int i;
2521
2522 shAssert(t_array != NULL);
2523 shAssert(numtriangles != NULL);
2524
2525 /* first, sort the array */
2526 sort_triangle_array(t_array, *numtriangles);
2527
2528 /*
2529 * now, find the first triangle with "ba" > AT_MATCH_RATIO and
2530 * re-set "numtriangles" to be just before it.
2531 *
2532 * if this would make "numtriangles" < 1, assert
2533 */
2534 for (i = (*numtriangles) - 1; i >= 0; i--) {
2535 if (t_array[i].ba <= AT_MATCH_RATIO) {
2536 break;
2537 }
2538 }
2539 *numtriangles = i;
2540 shAssert(*numtriangles >= 0);
2541 }
2542
2543
2544
2545 /************************************************************************
2546 *
2547 *
2548 * ROUTINE: make_vote_matrix
2549 *
2550 * DESCRIPTION:
2551 * Given two arrays of triangles, and the arrays of stars that make
2552 * up each set of triangles, try to match up triangles in the two
2553 * arrays. Triangles can be considered to match only when the
2554 * Euclidean distance in "triangle space", created from the two
2555 * coordinates "ba" and "ca", is within "max_radius". That is,
2556 * for two triangles to match, we must satisfy
2557 *
2558 * sqrt[ (t1.ba - t2.ba)^2 + (t1.ca - t2.ca)^2 ] <= max_radius
2559 *
2560 * Note that there may be more than one triangle from array A which
2561 * matches a particular triangle from array B! That's okay --
2562 * we treat any 2 which satisfy the above equation as "matched".
2563 * We rely upon the "vote_array" to weed out false matches.
2564 *
2565 * If "min_scale" and "max_scale" are not both -1, then disallow
2566 * any match for which the
2567 * ratio of triangles (indicated by "a_length" members)
2568 * is outside the given values.
2569 *
2570 * If "rotation_deg" and "tolerance_deg" are not both AT_MATCH_NOANGLE,
2571 * then disallow any match for which the two triangles
2572 * are not oriented at an angle of "rotation_deg" degrees
2573 * relative to each other (with a tolerance of "tolerance_deg" degrees).
2574 *
2575 * For each pair of triangles that matches, increment
2576 * the "vote" in each "vote cell" for each pair of matching
2577 * vertices.
2578 *
2579 * The "vote matrix" is a 2-D array of 'nbright'-by-'nbright'
2580 * integers. We allocate the array in this function, and
2581 * return a pointer to the array. Each cell in the array, vote[i][j],
2582 * contains the number of triangles in which
2583 *
2584 * star_array_A[i] matched star_array_B[j]
2585 *
2586 *
2587 * RETURN:
2588 * int ** pointer to new "vote matrix"
2589 *
2590 * </AUTO>
2591 */
2592
2593 static int **
make_vote_matrix(s_star * star_array_A,int num_stars_A,s_star * star_array_B,int num_stars_B,s_triangle * t_array_A,int num_triangles_A,s_triangle * t_array_B,int num_triangles_B,int nbright,double max_radius,double min_scale,double max_scale,double rotation_deg,double tolerance_deg)2594 make_vote_matrix
2595 (
2596 s_star *star_array_A, /* I: first array of stars */
2597 int num_stars_A, /* I: number of stars in star_array_A */
2598 s_star *star_array_B, /* I: second array of stars */
2599 int num_stars_B, /* I: number of stars in star_array_B */
2600 s_triangle *t_array_A, /* I: array of triangles from star_array_A */
2601 int num_triangles_A, /* I: number of triangles in t_array_A */
2602 s_triangle *t_array_B, /* I: array of triangles from star_array_B */
2603 int num_triangles_B, /* I: number of triangles in t_array_B */
2604 int nbright, /* I: consider at most this many stars */
2605 /* from each array; also the size */
2606 /* of the output "vote_matrix". */
2607 double max_radius, /* I: max radius in triangle-space allowed */
2608 /* for 2 triangles to be considered */
2609 /* a matching pair. */
2610 double min_scale, /* I: minimum permitted relative scale factor */
2611 /* if -1, any scale factor is allowed */
2612 double max_scale, /* I: maximum permitted relative scale factor */
2613 /* if -1, any scale factor is allowed */
2614 double rotation_deg, /* I: desired relative angle of coord systems (deg) */
2615 /* if AT_MATCH_NOANGLE, any orientation is allowed */
2616 double tolerance_deg /* I: allowed range of orientation angles (deg) */
2617 /* if AT_MATCH_NOANGLE, any orientation is allowed */
2618 )
2619 {
2620 int i, j, start_index;
2621 int **vote_matrix;
2622 double ba_A, ba_B, ca_A, ca_B, ba_min, ba_max;
2623 double rad2;
2624 double ratio;
2625 double actual_angle_deg;
2626 struct s_triangle *tri;
2627
2628 shAssert(star_array_A != NULL);
2629 shAssert(star_array_B != NULL);
2630 shAssert(t_array_A != NULL);
2631 shAssert(t_array_B != NULL);
2632 shAssert(nbright > 0);
2633 if (min_scale != -1) {
2634 shAssert((max_scale != -1) && (min_scale <= max_scale));
2635 }
2636 if (max_scale != -1) {
2637 shAssert((min_scale != -1) && (min_scale <= max_scale));
2638 }
2639
2640
2641 /* allocate and initialize the "vote_matrix" */
2642 vote_matrix = (int **) shMalloc(nbright*sizeof(int *));
2643 for (i = 0; i < nbright; i++) {
2644 vote_matrix[i] = (int *) shMalloc(nbright*sizeof(int));
2645 for (j = 0; j < nbright; j++) {
2646 vote_matrix[i][j] = 0;
2647 }
2648 }
2649
2650 /*
2651 * now, the triangles in "t_array_A" have been sorted by their "ba"
2652 * values. Therefore, we walk through the OTHER array, "t_array_B",
2653 * and for each triangle tri_B in it
2654 *
2655 * 1a. set ba_min = tri_B.ba - max_radius
2656 * 1b. set ba_max = tri_B.ba + max_radius
2657 *
2658 * We'll use these values to limit our selection from t_array_A.
2659 *
2660 * 2. find the first triangle in t_array_A which has
2661 * ba > ba_min
2662 * 3. starting there, step through t_array_A, calculating the
2663 * Euclidean distance between tri_B and
2664 * the current triangle from array A.
2665 * 4. stop stepping through t_array_A when we each a triangle
2666 * with ba > ba_max
2667 */
2668 rad2 = max_radius*max_radius;
2669 for (j = 0; j < num_triangles_B; j++) {
2670
2671 /*
2672 * make sure that this triangle doesn't have a vertex with index
2673 * greater than n_bright (because, if it did, we'd overwrite memory
2674 * when we tried to increment the vote_matrix array element).
2675 *
2676 * This is only a problem when called from "smallTrans", with
2677 * num_stars_A > nbright
2678 * or
2679 * num_stars_B > nbright
2680 */
2681 tri = &(t_array_B[j]);
2682 if ((tri->a_index >= nbright) || (tri->b_index >= nbright) ||
2683 (tri->c_index >= nbright)) {
2684 #ifdef DEBUG2
2685 printf("make_vote_matrix: skipping B triangle %d\n", j);
2686 #endif
2687 continue;
2688 }
2689
2690 #ifdef DEBUG2
2691 printf("make_vote_matrix: looking for matches to B %d\n", j);
2692 #endif
2693 ba_B = t_array_B[j].ba;
2694 ca_B = t_array_B[j].ca;
2695 ba_min = ba_B - max_radius;
2696 ba_max = ba_B + max_radius;
2697 #ifdef DEBUG2
2698 printf(" ba_min = %7.3f ba_max = %7.3f\n", ba_min, ba_max);
2699 #endif
2700
2701 start_index = find_ba_triangle(t_array_A, num_triangles_A, ba_min);
2702 for (i = start_index; i < num_triangles_A; i++) {
2703
2704 /*
2705 * again, skip any triangle which has a vertex with ID > nbright
2706 */
2707 tri = &(t_array_A[i]);
2708 if ((tri->a_index >= nbright) || (tri->b_index >= nbright) ||
2709 (tri->c_index >= nbright)) {
2710 #ifdef DEBUG2
2711 printf("make_vote_matrix: skipping A triangle %d\n", i);
2712 #endif
2713 continue;
2714 }
2715
2716 #ifdef DEBUG2
2717 printf(" looking at A %d\n", i);
2718 #endif
2719 ba_A = t_array_A[i].ba;
2720 ca_A = t_array_A[i].ca;
2721
2722 /* check to see if we can stop looking through A yet */
2723 if (ba_A > ba_max) {
2724 break;
2725 }
2726
2727 if ((ba_A - ba_B)*(ba_A - ba_B) + (ca_A - ca_B)*(ca_A - ca_B) < rad2) {
2728
2729 /*
2730 * check the ratio of lengths of side "a", and discard this
2731 * candidate if its outside the allowed range
2732 */
2733 if (min_scale != -1) {
2734 ratio = t_array_A[i].a_length/t_array_B[j].a_length;
2735 if (ratio < min_scale || ratio > max_scale) {
2736 continue;
2737 }
2738 }
2739
2740 /*
2741 * check the relative orientations of the triangles.
2742 * If they don't match the desired rotation angle,
2743 * discard this match.
2744 */
2745 if (rotation_deg != AT_MATCH_NOANGLE) {
2746 if (is_desired_rotation(&(t_array_A[i]), &(t_array_B[j]),
2747 rotation_deg, tolerance_deg, &actual_angle_deg) == 0) {
2748 continue;
2749 }
2750 }
2751
2752
2753 /* we have a (possible) match! */
2754 #ifdef DEBUG2
2755 ratio = t_array_A[i].a_length/t_array_B[j].a_length;
2756 printf(" match! A: (%6.3f, %6.3f) B: (%6.3f, %6.3f) ratio %9.4e angle %9.4e\n",
2757 ba_A, ca_A, ba_B, ca_B, ratio, actual_angle_deg);
2758 #endif
2759 /*
2760 * increment the vote_matrix cell for each matching pair
2761 * of stars, one at each vertex
2762 */
2763 vote_matrix[t_array_A[i].a_index][t_array_B[j].a_index]++;
2764 vote_matrix[t_array_A[i].b_index][t_array_B[j].b_index]++;
2765 vote_matrix[t_array_A[i].c_index][t_array_B[j].c_index]++;
2766
2767 }
2768 }
2769 }
2770
2771 #ifdef DEBUG
2772 print_vote_matrix(vote_matrix, nbright);
2773 #endif
2774
2775
2776 return(vote_matrix);
2777 }
2778
2779
2780 /************************************************************************
2781 *
2782 *
2783 * ROUTINE: print_vote_matrix
2784 *
2785 * DESCRIPTION:
2786 * Print out the "vote_matrix" in a nice format.
2787 *
2788 * For debugging purposes.
2789 *
2790 * RETURN:
2791 * nothing
2792 *
2793 * </AUTO>
2794 */
2795
2796 #ifdef DEBUG
2797
2798 static void
print_vote_matrix(int ** vote_matrix,int numcells)2799 print_vote_matrix
2800 (
2801 int **vote_matrix, /* I: the 2-D array we'll print out */
2802 int numcells /* I: number of cells in each row and col of matrix */
2803 )
2804 {
2805 int i, j;
2806
2807 printf("here comes vote matrix\n");
2808 for (i = 0; i < numcells; i++) {
2809 for (j = 0; j < numcells; j++) {
2810 printf(" %3d", vote_matrix[i][j]);
2811 }
2812 printf("\n");
2813 }
2814 }
2815
2816 #endif /* DEBUG */
2817
2818
2819 /************************************************************************
2820 *
2821 *
2822 * ROUTINE: top_vote_getters
2823 *
2824 * DESCRIPTION:
2825 * Given a vote_matrix which has been filled in,
2826 * which has 'num' rows and columns, we need to pick the
2827 * top 'num' vote-getters. We call 'top_vote_getters'
2828 * and are given, in its output arguments, pointers to three
2829 * arrays, each of which has 'num' elements pertaining
2830 * to a matched pair of STARS:
2831 *
2832 * winner_votes[] number of votes of winners, in descending order
2833 * winner_index_A[] index of star in star_array_A
2834 * winner_index_B[] index of star in star_array_B
2835 *
2836 * Thus, the pair of stars which matched in the largest number
2837 * of triangles will be
2838 *
2839 * star_array_A[winner_index_A[0]] from array A
2840 * star_array_B[winner_index_A[0]] from array B
2841 *
2842 * and the pair of stars which matched in the second-largest number
2843 * of triangles will be
2844 *
2845 * star_array_A[winner_index_A[1]] from array A
2846 * star_array_B[winner_index_A[1]] from array B
2847 *
2848 * and so on.
2849 *
2850 * RETURN:
2851 * SH_SUCCESS if all goes well
2852 * SH_GENERIC_ERROR if not
2853 *
2854 * </AUTO>
2855 */
2856
2857 static int
top_vote_getters(int ** vote_matrix,int num,int ** winner_votes,int ** winner_index_A,int ** winner_index_B)2858 top_vote_getters
2859 (
2860 int **vote_matrix, /* I: the 2-D array, already filled in */
2861 int num, /* I: # of rows and cols in vote_matrix */
2862 /* also the number of elements in the next */
2863 /* three output arrays */
2864 int **winner_votes, /* O: create this array of # of votes for the */
2865 /* 'num' cells with the most votes */
2866 int **winner_index_A, /* O: create this array of index into star array A */
2867 /* of the 'num' cells with most votes */
2868 int **winner_index_B /* O: create this array of index into star array B */
2869 /* of the 'num' cells with most votes */
2870 )
2871 {
2872 int i, j, k, l;
2873 int *w_votes; /* local ptr to (*winner_votes), for convenience */
2874 int *w_index_A; /* local ptr to (*winner_index_A), for convenience */
2875 int *w_index_B; /* local ptr to (*winner_index_B), for convenience */
2876
2877 /* first, create the output arrays */
2878 *winner_votes = (int *) shMalloc(num*sizeof(int));
2879 *winner_index_A = (int *) shMalloc(num*sizeof(int));
2880 *winner_index_B = (int *) shMalloc(num*sizeof(int));
2881
2882 /* this will simplify code inside this function */
2883 w_votes = *winner_votes;
2884 w_index_A = *winner_index_A;
2885 w_index_B = *winner_index_B;
2886
2887 /*
2888 * initialize all elements of the output arrays. Use -1 as the
2889 * index in "w_index" arrays, to indicate an empty place
2890 * with no real winner.
2891 */
2892 for (i = 0; i < num; i++) {
2893 w_votes[i] = 0;
2894 w_index_A[i] = -1;
2895 w_index_B[i] = -1;
2896 }
2897
2898 /*
2899 * now walk through the vote_matrix, using insertion sort to place
2900 * a cell into the "winner" arrays if it has more votes than the
2901 * least popular winner so far (i.e. w_votes[num - 1])
2902 */
2903 for (i = 0; i < num; i++) {
2904 for (j = 0; j < num; j++) {
2905 if (vote_matrix[i][j] > w_votes[num - 1]) {
2906
2907 /* have to insert this cell's values into the winner arrays */
2908 for (k = 0; k < num; k++) {
2909 if (vote_matrix[i][j] > w_votes[k]) {
2910
2911 /* move all other winners down one place */
2912 for (l = num - 2; l >= k; l--) {
2913 w_votes[l + 1] = w_votes[l];
2914 w_index_A[l + 1] = w_index_A[l];
2915 w_index_B[l + 1] = w_index_B[l];
2916 }
2917 /* insert the new item in its place */
2918 w_votes[k] = vote_matrix[i][j];
2919 w_index_A[k] = i;
2920 w_index_B[k] = j;
2921 break;
2922 }
2923 }
2924 }
2925 }
2926 }
2927
2928 #ifdef DEBUG
2929 printf(" in top_vote_getters, we have top %d \n", num);
2930 for (i = 0; i < num; i++) {
2931 printf(" index_A %4d index_B %4d votes %4d\n",
2932 w_index_A[i], w_index_B[i], w_votes[i]);
2933 }
2934 #endif
2935
2936 return(SH_SUCCESS);
2937 }
2938
2939
2940
2941 /************************************************************************
2942 *
2943 *
2944 * ROUTINE: calc_trans
2945 *
2946 * DESCRIPTION:
2947 * Given a set of "nbright" matched pairs of stars, which we can
2948 * extract from the "winner_index" and "star_array" arrays,
2949 * figure out a TRANS structure which takes coordinates of
2950 * objects in set A and transforms then into coords for set B.
2951 * A TRANS contains 6, 12, or 16 coefficients in equations like this:
2952 *
2953 * if linear terms only:
2954 *
2955 * x' = A + B*x + C*y
2956 * y' = D + E*x + F*y
2957 *
2958 * if linear plus quadratic terms,
2959 *
2960 * x' = A + Bx + Cy + Dxx + Exy + Fyy
2961 * y' = G + Hx + Iy + Jxx + Kxy + Lyy
2962 *
2963 * if linear plus quadratic plus cubic,
2964 *
2965 * x' = A + Bx + Cy + Dxx + Exy + Fyy + Gx(xx+yy) + Hy(xx+yy)
2966 * y' = I + Jx + Ky + Lxx + Mxy + Nyy + Ox(xx+yy) + Py(xx+yy)
2967 *
2968 * where (x,y) are coords in set A and (x',y') are corresponding
2969 * coords in set B.
2970 *
2971 * This function simply checks the value of the TRANS 'order' field,
2972 * and calls the appropriate function to do the actual work.
2973 *
2974 *
2975 * RETURN:
2976 * SH_SUCCESS if all goes well
2977 * SH_GENERIC_ERROR if we can't find a solution
2978 *
2979 * </AUTO>
2980 */
2981
2982 static int
calc_trans(int nbright,s_star * star_array_A,int num_stars_A,s_star * star_array_B,int num_stars_B,int * winner_votes,int * winner_index_A,int * winner_index_B,TRANS * trans)2983 calc_trans
2984 (
2985 int nbright, /* I: max number of stars we use in calculating */
2986 /* the transformation; we may cut down to */
2987 /* a more well-behaved subset. */
2988 s_star *star_array_A, /* I: first array of s_star structure we match */
2989 /* the output TRANS takes their coords */
2990 /* into those of array B */
2991 int num_stars_A, /* I: total number of stars in star_array_A */
2992 s_star *star_array_B, /* I: second array of s_star structure we match */
2993 int num_stars_B, /* I: total number of stars in star_array_B */
2994 int *winner_votes, /* I: number of votes gotten by the top 'nbright' */
2995 /* matched pairs of stars */
2996 int *winner_index_A, /* I: index into "star_array_A" of top */
2997 /* vote-getters */
2998 int *winner_index_B, /* I: index into "star_array_B" of top */
2999 /* vote-getters */
3000 TRANS *trans /* O: place solved coefficients into this */
3001 /* existing structure's fields */
3002 )
3003 {
3004
3005 /*
3006 * using the trans->order value, call the appropriate function
3007 */
3008 switch (trans->order) {
3009 case AT_TRANS_LINEAR:
3010 if (calc_trans_linear(nbright,
3011 star_array_A, num_stars_A,
3012 star_array_B, num_stars_B,
3013 winner_votes, winner_index_A, winner_index_B,
3014 trans) != SH_SUCCESS) {
3015 shError("calc_trans: calc_trans_linear returns with error");
3016 return(SH_GENERIC_ERROR);
3017 }
3018 break;
3019
3020 case AT_TRANS_QUADRATIC:
3021 if (calc_trans_quadratic(nbright,
3022 star_array_A, num_stars_A,
3023 star_array_B, num_stars_B,
3024 winner_votes, winner_index_A, winner_index_B,
3025 trans) != SH_SUCCESS) {
3026 shError("calc_trans: calc_trans_quadratic returns with error");
3027 return(SH_GENERIC_ERROR);
3028 }
3029 break;
3030
3031 case AT_TRANS_CUBIC:
3032 if (calc_trans_cubic(nbright,
3033 star_array_A, num_stars_A,
3034 star_array_B, num_stars_B,
3035 winner_votes, winner_index_A, winner_index_B,
3036 trans) != SH_SUCCESS) {
3037 shError("calc_trans: calc_trans_cubic returns with error");
3038 return(SH_GENERIC_ERROR);
3039 }
3040 break;
3041
3042 default:
3043 shFatal("calc_trans: called with invalid trans->order %d \n",
3044 trans->order);
3045 break;
3046 }
3047
3048 return(SH_SUCCESS);
3049 }
3050
3051
3052
3053 /************************************************************************
3054 *
3055 *
3056 * ROUTINE: alloc_matrix
3057 *
3058 * DESCRIPTION:
3059 * Allocate space for an NxN matrix of double values,
3060 * return a pointer to the new matrix.
3061 *
3062 * RETURNS:
3063 * double ** pointer to new matrix
3064 *
3065 *
3066 * </AUTO>
3067 */
3068
3069 static double **
alloc_matrix(int n)3070 alloc_matrix
3071 (
3072 int n /* I: number of elements in each row and col */
3073 )
3074 {
3075 int i;
3076 double **matrix;
3077
3078 matrix = (double **) shMalloc(n*sizeof(double *));
3079 for (i = 0; i < n; i++) {
3080 matrix[i] = (double *) shMalloc(n*sizeof(double));
3081 }
3082
3083 return(matrix);
3084 }
3085
3086
3087 /************************************************************************
3088 *
3089 *
3090 * ROUTINE: free_matrix
3091 *
3092 * DESCRIPTION:
3093 * Free the space allocated for the given nxn matrix.
3094 *
3095 * RETURNS:
3096 * nothing
3097 *
3098 * </AUTO>
3099 */
3100
3101 static void
free_matrix(double ** matrix,int n)3102 free_matrix
3103 (
3104 double **matrix, /* I: pointer to 2-D array to be freed */
3105 int n /* I: number of elements in each row and col */
3106 )
3107 {
3108 int i;
3109
3110 for (i = 0; i < n; i++) {
3111 shFree(matrix[i]);
3112 }
3113 shFree(matrix);
3114 }
3115
3116
3117 /************************************************************************
3118 *
3119 *
3120 * ROUTINE: print_matrix
3121 *
3122 * DESCRIPTION:
3123 * print out a nice picture of the given matrix.
3124 *
3125 * For debugging purposes.
3126 *
3127 * RETURNS:
3128 * nothing
3129 *
3130 * </AUTO>
3131 */
3132
3133 #ifdef DEBUG
3134
3135 static void
print_matrix(double ** matrix,int n)3136 print_matrix
3137 (
3138 double **matrix, /* I: pointer to 2-D array to be printed */
3139 int n /* I: number of elements in each row and col */
3140 )
3141 {
3142 int i, j;
3143
3144 for (i = 0; i < n; i++) {
3145 for (j = 0; j < n; j++) {
3146 printf(" %12.5e", matrix[i][j]);
3147 }
3148 printf("\n");
3149 }
3150 }
3151
3152 #endif /* DEBUG */
3153
3154
3155
3156
3157 /*
3158 * check to see if my versions of NR routines have bugs.
3159 * Try to invert a matrix.
3160 *
3161 * debugging only.
3162 */
3163
3164 #ifdef DEBUG3
3165
3166 static void
test_routine(void)3167 test_routine (void)
3168 {
3169 int i, j, k, n;
3170 int *permutations;
3171 double **matrix1, **matrix2, **inverse;
3172 double *vector;
3173 double *col;
3174 double sum;
3175
3176 fflush(stdout);
3177 fflush(stderr);
3178 n = 2;
3179 matrix1 = (double **) shMalloc(n*sizeof(double *));
3180 matrix2 = (double **) shMalloc(n*sizeof(double *));
3181 inverse = (double **) shMalloc(n*sizeof(double *));
3182 vector = (double *) shMalloc(n*sizeof(double));
3183 for (i = 0; i < n; i++) {
3184 matrix1[i] = (double *) shMalloc(n*sizeof(double));
3185 matrix2[i] = (double *) shMalloc(n*sizeof(double));
3186 inverse[i] = (double *) shMalloc(n*sizeof(double));
3187 }
3188 permutations = (int *) shMalloc(n*sizeof(int));
3189 col = (double *) shMalloc(n*sizeof(double));
3190
3191
3192 /* fill the matrix */
3193 matrix1[0][0] = 1.0;
3194 matrix1[0][1] = 2.0;
3195 matrix1[1][0] = 3.0;
3196 matrix1[1][1] = 4.0;
3197
3198 /* fill the vector */
3199 for (i = 0; i < n; i++) {
3200 vector[i] = 0;
3201 }
3202
3203 /* copy matrix1 into matrix2, so we can compare them later */
3204 for (i = 0; i < n; i++) {
3205 for (j = 0; j < n; j++) {
3206 matrix2[i][j] = matrix1[i][j];
3207 }
3208 }
3209
3210 /* now check */
3211 printf(" here comes original matrix \n");
3212 print_matrix(matrix1, n);
3213
3214
3215 /* now invert matrix1 */
3216 for (i = 0; i < n; i++) {
3217 for (j = 0; j < n; j++) {
3218 inverse[i][j] = matrix1[i][j];
3219 }
3220 }
3221 gauss_matrix(inverse, n, vector);
3222
3223 /* now check */
3224 printf(" here comes inverse matrix \n");
3225 print_matrix(inverse, n);
3226
3227 /* find out if the product of "inverse" and "matrix2" is identity */
3228 sum = 0.0;
3229 for (i = 0; i < n; i++) {
3230 for (j = 0; j < n; j++) {
3231 for (k = 0; k < n; k++) {
3232 sum += inverse[i][k]*matrix2[k][j];
3233 }
3234 matrix1[i][j] = sum;
3235 sum = 0.0;
3236 }
3237 }
3238
3239 printf(" here comes what we hope is identity matrix \n");
3240 print_matrix(matrix1, n);
3241
3242 fflush(stdout);
3243 fflush(stderr);
3244 }
3245
3246 #endif /* DEBUG3 */
3247
3248
3249
3250
3251
3252
3253
3254
3255 /************************************************************************
3256 *
3257 *
3258 * ROUTINE: iter_trans
3259 *
3260 * DESCRIPTION:
3261 * We want to find a TRANS structures that takes coords of objects in
3262 * set A and transforms to coords of objects in set B. We have a
3263 * a subset of 'nmatched' candidates for matched pairs of points.
3264 * However, some of these may be false matches. Here's how we try
3265 * to eliminate them, and use all remaining true matches to derive
3266 * the transformation.
3267 *
3268 * 1. start with nbright matched candidate pairs of points
3269 * 2. choose N best pairs
3270 * if recalc_flag == RECALC_NO, set N = AT_MATCH_STARTN
3271 * if recalc_flag == RECALC_YES, set N = nbright
3272 * (assert that N >= AT_MATCH_REQUIRE)
3273 * 3. set Nr = N
3274 * 4. calculate a TRANS structure using the best Nr points
3275 * (where "best" means "highest in winner_index" arrays)
3276 * 5. transform all Nr points from coords in A to coords in B
3277 * 6. calculate Euclidean square-of-distance between all Nr points
3278 * in coord system B
3279 * 7. sort these Euclidean values
3280 * 8. pick the AT_MATCH_PERCENTILE'th value from the sorted array
3281 * (call it "sigma")
3282 * 9. let Nb = number of candidate matched pairs which have
3283 * square-of-distance > AT_MATCH_NSIGMA*sigma
3284 * 10. if Nb == 0, we're done -- quit
3285 * 11. if Nb > 0,
3286 * remove all Nb candidates from matched pair arrays
3287 * set Nr = Nr - Nb
3288 * go to step 4
3289 *
3290 * Note that if we run out of candidate pairs, so that Nr < AT_MATCH_REQUIRE,
3291 * we print an error message and return SH_GENERIC_ERROR.
3292 *
3293 * The "recalc_flag" is used to distinguish two cases:
3294 * if RECALC_NO, then we are calling 'iter_trans()' with a bunch of
3295 * matches which probably contain some bad ones.
3296 * In order to prevent the bad ones from ruining the
3297 * initial calculation, we pick only the few best
3298 * on the first iteration.
3299 * if RECALC_YES, we are calling 'iter_trans()' with a set of matches
3300 * which have already passed a test: namely, they are
3301 * based on a previously-determined TRANS and all matches
3302 * are within 'matchrad' in the coord system of list B.
3303 * In this case, we start out using all the matched
3304 * pairs in the very first iteration.
3305 *
3306 *
3307 * RETURNS:
3308 * SH_SUCCESS if we were able to determine a good TRANS
3309 * SH_GENERIC_ERROR if we couldn't
3310 *
3311 * </AUTO>
3312 */
3313
3314 static int
iter_trans(int nbright,s_star * star_array_A,int num_stars_A,s_star * star_array_B,int num_stars_B,int * winner_votes,int * winner_index_A,int * winner_index_B,int recalc_flag,int max_iterations,double halt_sigma,TRANS * trans)3315 iter_trans
3316 (
3317 int nbright, /* I: max number of stars we use in calculating */
3318 /* the transformation; we may cut down to */
3319 /* a more well-behaved subset. */
3320 s_star *star_array_A, /* I: first array of s_star structure we match */
3321 /* the output TRANS takes their coords */
3322 /* into those of array B */
3323 int num_stars_A, /* I: total number of stars in star_array_A */
3324 s_star *star_array_B, /* I: second array of s_star structure we match */
3325 int num_stars_B, /* I: total number of stars in star_array_B */
3326 int *winner_votes, /* I: number of votes gotten by the top 'nbright' */
3327 /* matched pairs of stars */
3328 /* We may modify this array */
3329 int *winner_index_A, /* I: index into "star_array_A" of top */
3330 /* vote-getters */
3331 /* We may modify this array */
3332 int *winner_index_B, /* I: index into "star_array_B" of top */
3333 /* vote-getters */
3334 /* We may modify this array */
3335 int recalc_flag, /* I: should we use only a few best pairs for */
3336 /* the first iteration, or all? */
3337 int max_iterations, /* I: iterate at most this many times. If we
3338 /* reach this limit, stop iterating */
3339 /* and declare success */
3340 double halt_sigma, /* I: if the residuals from solution drop to */
3341 /* this level, stop iterating and */
3342 /* declare success */
3343 TRANS *trans /* O: place solved coefficients into this */
3344 /* existing structure's fields */
3345 )
3346 {
3347 int i, j;
3348 int nr; /* number of matched pairs remaining in solution */
3349 int nb; /* number of bad pairs in any iteration */
3350 int initial_pairs;
3351 int is_ok;
3352 int required_pairs, start_pairs;
3353 int iters_so_far;
3354 double *dist2, *dist2_sorted;
3355 double xdiff, ydiff;
3356 double sigma;
3357 double max_dist2;
3358 double newx, newy;
3359 s_star *sa, *sb;
3360 s_star *a_prime; /* will hold transformed version of stars in set A */
3361
3362
3363 /*
3364 * set some variables depending on the order of the fit to be
3365 * performed.
3366 */
3367 switch (trans->order) {
3368 case AT_TRANS_LINEAR:
3369 required_pairs = AT_MATCH_REQUIRE_LINEAR;
3370 start_pairs = AT_MATCH_STARTN_LINEAR;
3371 break;
3372 case AT_TRANS_QUADRATIC:
3373 required_pairs = AT_MATCH_REQUIRE_QUADRATIC;
3374 start_pairs = AT_MATCH_STARTN_QUADRATIC;
3375 break;
3376 case AT_TRANS_CUBIC:
3377 required_pairs = AT_MATCH_REQUIRE_CUBIC;
3378 start_pairs = AT_MATCH_STARTN_CUBIC;
3379 break;
3380 default:
3381 shFatal("iter_trans: invalid trans->order %d \n", trans->order);
3382 break;
3383 }
3384
3385
3386 if (nbright < required_pairs) {
3387 #ifdef DEBUG
3388 printf("iter_trans: only %d items supplied, need %d\n",
3389 nbright, required_pairs);
3390 #endif
3391 return(SH_GENERIC_ERROR);
3392 }
3393
3394
3395 shAssert(star_array_A != NULL);
3396 shAssert(star_array_B != NULL);
3397 shAssert(winner_votes != NULL);
3398 shAssert(winner_index_A != NULL);
3399 shAssert(winner_index_B != NULL);
3400 shAssert(trans != NULL);
3401
3402 /* these should already have been checked, but it doesn't hurt */
3403 shAssert(num_stars_A >= nbright);
3404 shAssert(num_stars_A >= nbright);
3405
3406 /*
3407 * make a first guess at TRANS;
3408 * use all the pairs, if we are sure they are "safe",
3409 * or only the best 'start_pairs', if we're not sure
3410 */
3411 if (recalc_flag == RECALC_YES) {
3412 initial_pairs = nbright;
3413 } else {
3414 initial_pairs = start_pairs;
3415 }
3416 #ifdef DEBUG
3417 printf(" on initial calc, use %d pairs\n", initial_pairs);
3418 #endif
3419 if (calc_trans(initial_pairs,
3420 star_array_A, num_stars_A,
3421 star_array_B, num_stars_B,
3422 winner_votes, winner_index_A, winner_index_B,
3423 trans) != SH_SUCCESS) {
3424 shError("iter_trans: calc_trans returns with error");
3425 return(SH_GENERIC_ERROR);
3426 }
3427 #ifdef DEBUG
3428 printf(" here comes initial TRANS\n");
3429 print_trans(trans);
3430 #endif
3431
3432 /*
3433 * Now, we are going to enter the iteration with a set of the "best"
3434 * matched pairs. Recall that
3435 * "winner_index" arrays are already sorted in decreasing order
3436 * of goodness, so that "winner_index_A[0]" is the best.
3437 * As we iterate, we may discard some matches, and then 'nr' will
3438 * get smaller. It must always be more than AT_MATCH_REQUIRE,
3439 * or else 'calc_trans' will fail.
3440 */
3441 nr = nbright;
3442
3443 /*
3444 * We're going to need an array of (at most) 'nbright' stars
3445 * which hold the coordinates of stars in set A, after they've
3446 * been transformed into coordinates system of set B.
3447 */
3448 a_prime = (s_star *) shMalloc(nbright*sizeof(s_star));
3449
3450 /*
3451 * And this will be an array to hold the Euclidean square-of-distance
3452 * between a transformed star from set A and its partner from set B.
3453 *
3454 * "dist2_sorted" is a copy of the array which we'll sort .. but we need
3455 * to keep the original order, too.
3456 */
3457 dist2 = (double *) shMalloc(nbright*sizeof(double));
3458 dist2_sorted = (double *) shMalloc(nbright*sizeof(double));
3459
3460 /*
3461 * we don't allow any candidate matches which cause the stars to
3462 * differ by more than this much in the common coord system.
3463 */
3464 max_dist2 = AT_MATCH_MAXDIST*AT_MATCH_MAXDIST;
3465
3466 /*
3467 * now, we enter a loop that may execute several times.
3468 * We calculate the transformation for current 'nr' best points,
3469 * then check to see if we should throw out any matches because
3470 * the resulting transformed coordinates are too discrepant.
3471 * We break out of this loop near the bottom, with a status
3472 * provided by "is_ok"
3473 *
3474 * is_ok = 1 all went well, can return success
3475 * is_ok = 0 we failed for some reason.
3476 */
3477 is_ok = 1;
3478
3479 iters_so_far = 0;
3480 while (iters_so_far < max_iterations) {
3481
3482 #ifdef DEBUG
3483 printf("iter_trans: at top of loop, nr=%4d iters_so_far=%4d\n",
3484 nr, iters_so_far);
3485 #endif
3486
3487 nb = 0;
3488
3489 /*
3490 * apply the TRANS to the A stars in all 'nr' matched pairs.
3491 * we make a new set of s_stars with the transformed coordinates,
3492 * called "a_prime".
3493 */
3494 for (i = 0; i < nr; i++) {
3495 sa = &(star_array_A[winner_index_A[i]]);
3496 if (calc_trans_coords(sa, trans, &newx, &newy) != SH_SUCCESS) {
3497 shError("iter_trans: calc_trans_coords fails");
3498 return(SH_GENERIC_ERROR);
3499 }
3500 a_prime[i].x = newx;
3501 a_prime[i].y = newy;
3502 }
3503
3504
3505 /*
3506 * calculate the square-of-distance between a transformed star
3507 * (from set A) and its partner from set B, in the coordinate system
3508 * of set B.
3509 */
3510 for (i = 0; i < nr; i++) {
3511 sb = &(star_array_B[winner_index_B[i]]);
3512 xdiff = a_prime[i].x - sb->x;
3513 ydiff = a_prime[i].y - sb->y;
3514 dist2[i] = (xdiff*xdiff + ydiff*ydiff);
3515 dist2_sorted[i] = dist2[i];
3516 #ifdef DEBUG
3517 printf(" match %3d (%12.5e,%12.5e) vs. (%12.5e,%12.5e) d2=%12.6e\n",
3518 i, a_prime[i].x, a_prime[i].y, sb->x, sb->y, dist2[i]);
3519 #endif
3520 }
3521
3522 /*
3523 * sort the array of square-of-distances
3524 */
3525 qsort((char *) dist2_sorted, nr, sizeof(double), (PFI) compare_double);
3526
3527
3528 /*
3529 * now, check to see if any matches have dist2 > max_dist2.
3530 * If so,
3531 *
3532 * - remove them from the winner_votes and winner_index arrays
3533 * - decrement 'nr'
3534 * - also decrement the loop counter 'i', because we're going
3535 * to move up all items in the "winner" and "dist2" arrays
3536 * as we discard the bad match
3537 * - increment 'nb'
3538 */
3539 for (i = 0; i < nr; i++) {
3540 if (dist2[i] > max_dist2) {
3541
3542 /*
3543 * remove the entry for the "bad" match from the "winner" arrays
3544 * and from the "dist2" array
3545 */
3546 #ifdef DEBUG
3547 printf(" removing old match with d2=%9.4e\n", dist2[i]);
3548 #endif
3549 for (j = i + 1; j < nr; j++) {
3550 winner_votes[j - 1] = winner_votes[j];
3551 winner_index_A[j - 1] = winner_index_A[j];
3552 winner_index_B[j - 1] = winner_index_B[j];
3553 dist2[j - 1] = dist2[j];
3554 }
3555
3556 /*
3557 * and modify our counters of "remaining good matches" and
3558 * "bad matches this time", too.
3559 */
3560 nr--; /* one fewer good match remains */
3561 nb++; /* one more bad match during this iteration */
3562
3563 /*
3564 * and decrement 'i', too, since we must moved element
3565 * i+1 to the place i used to be, and we must check _it_.
3566 */
3567 i--;
3568 }
3569 }
3570 #ifdef DEBUG
3571 printf(" nr now %4d, nb now %4d\n", nr, nb);
3572 #endif
3573
3574
3575 /*
3576 * pick the square-of-distance which occurs at the AT_MATCH_PERCENTILE
3577 * place in the sorted array. Call this value "sigma". We'll clip
3578 * any matches that are more than AT_MATCH_NSIGMA*"sigma".
3579 *
3580 * However, if we have fewer than 2 objects, don't bother with this
3581 * step -- just set "sigma" equal to 0 and prepare for later
3582 * failure....
3583 */
3584 if (nr < 2) {
3585 sigma = 0.0;
3586 #ifdef DEBUG
3587 printf(" sigma = %10.5e (only %d matches) \n", sigma, nr);
3588 #endif
3589 }
3590 else {
3591 sigma = find_percentile(dist2_sorted, nr, (double) AT_MATCH_PERCENTILE);
3592 #ifdef DEBUG
3593 printf(" sigma = %10.5e\n", sigma);
3594 #endif
3595 }
3596
3597 /*
3598 * If the current "sigma" value is less than the "halt_sigma" value,
3599 * then we have succeeded. Stop iterating.
3600 */
3601 if (sigma <= halt_sigma) {
3602 #ifdef DEBUG
3603 printf(" SUCCESS sigma = %10.5e < halt_sigma %10.5e \n",
3604 sigma, halt_sigma);
3605 #endif
3606 is_ok = 1;
3607 break;
3608 }
3609
3610 /*
3611 * now, check to see if any matches have dist2 > AT_MATCH_NSIGMA*sigma.
3612 * If so,
3613 *
3614 * - remove them from the winner_votes and winner_index arrays
3615 * - decrement 'nr'
3616 * - also decrement the loop counter 'i', because we're going
3617 * to move up all items in the "winner" and "dist2" arrays
3618 * as we discard the bad match
3619 * - increment 'nb'
3620 */
3621 for (i = 0; i < nr; i++) {
3622 if (dist2[i] > AT_MATCH_NSIGMA*sigma) {
3623
3624 /*
3625 * remove the entry for the "bad" match from the "winner" arrays
3626 * and from the "dist2" array
3627 */
3628 #ifdef DEBUG
3629 printf(" removing old match with d2=%9.4e\n", dist2[i]);
3630 #endif
3631 for (j = i + 1; j < nr; j++) {
3632 winner_votes[j - 1] = winner_votes[j];
3633 winner_index_A[j - 1] = winner_index_A[j];
3634 winner_index_B[j - 1] = winner_index_B[j];
3635 dist2[j - 1] = dist2[j];
3636 }
3637
3638 /*
3639 * and modify our counters of "remaining good matches" and
3640 * "bad matches this time", too.
3641 */
3642 nr--; /* one fewer good match remains */
3643 nb++; /* one more bad match during this iteration */
3644
3645 /*
3646 * and decrement 'i', too, since we must moved element
3647 * i+1 to the place i used to be, and we must check _it_.
3648 */
3649 i--;
3650 }
3651 }
3652 #ifdef DEBUG
3653 printf(" nr now %4d, nb now %4d\n", nr, nb);
3654 #endif
3655
3656
3657 /*
3658 * Okay, let's evaluate what has happened so far:
3659 * - if nb == 0, then all remaining matches are good
3660 * - if nb > 0, we need to iterate again
3661 * - if nr < required_pairs, we've thrown out too many points,
3662 * and must quit in shame
3663 */
3664 if (nb == 0) {
3665 #ifdef DEBUG
3666 printf(" SUCCESS nb = 0, no more pairs to discard \n");
3667 #endif
3668 is_ok = 1;
3669 break;
3670 }
3671
3672 if (nr < required_pairs) {
3673 shDebug(AT_MATCH_ERRLEVEL,
3674 "iter_trans: only %d points remain, fewer than %d required",
3675 nr, required_pairs);
3676 is_ok = 0;
3677 break;
3678 }
3679
3680
3681 /*
3682 * calculate the TRANS for the remaining set of matches
3683 */
3684 #ifdef DEBUG
3685 printf(" on this iter, use %d pairs\n", nr);
3686 #endif
3687 if (calc_trans(nr, star_array_A, num_stars_A,
3688 star_array_B, num_stars_B,
3689 winner_votes, winner_index_A, winner_index_B,
3690 trans) != SH_SUCCESS) {
3691 shError("iter_trans: calc_trans returns with error");
3692 return(SH_GENERIC_ERROR);
3693 }
3694
3695 #ifdef DEBUG
3696 printf(" here comes latest TRANS\n");
3697 print_trans(trans);
3698 #endif
3699
3700 iters_so_far++;
3701
3702 }
3703
3704 if (iters_so_far == max_iterations) {
3705 #ifdef DEBUG
3706 printf(" SUCCESS(?): iters_so_far %d = max_iterations\n",
3707 iters_so_far);
3708 #endif
3709 }
3710
3711 /*
3712 * Here we summarize the result of our work in two of the
3713 * elements of the TRANS structure:
3714 * trans->nr = number of pairs used to find transformation
3715 * trans->sig = stdev of separation of matching pairs,
3716 * in units of coord system B
3717 */
3718 trans->nr = nr;
3719 trans->sig = find_percentile(dist2_sorted, nr, ONE_STDEV_PERCENTILE);
3720
3721 /*
3722 * free up the arrays we allocated
3723 */
3724 shFree(a_prime);
3725 shFree(dist2);
3726 shFree(dist2_sorted);
3727
3728 /*
3729 * and decide whether we succeeded, or failed
3730 */
3731 if (is_ok == 0) {
3732 return(SH_GENERIC_ERROR);
3733 }
3734 else {
3735 return(SH_SUCCESS);
3736 }
3737 }
3738
3739
3740
3741 /************************************************************************
3742 *
3743 *
3744 * ROUTINE: compare_double
3745 *
3746 * DESCRIPTION:
3747 * Given pointers to two double numbers, return the comparison.
3748 * Used by "iter_trans"
3749 *
3750 * RETURN:
3751 * 1 if first double is larger than second
3752 * 0 if the two are equal
3753 * -1 if first double is smaller than second
3754 *
3755 * </AUTO>
3756 */
3757
3758 static int
compare_double(double * f1,double * f2)3759 compare_double
3760 (
3761 double *f1, /* I: compare size of FIRST double value */
3762 double *f2 /* I: ... with SECOND double value */
3763 )
3764 {
3765 shAssert((f1 != NULL) && (f2 != NULL));
3766
3767 if (*f1 > *f2) {
3768 return(1);
3769 }
3770 if (*f1 < *f2) {
3771 return(-1);
3772 }
3773 return(0);
3774 }
3775
3776
3777 /************************************************************************
3778 *
3779 *
3780 * ROUTINE: find_percentile
3781 *
3782 * DESCRIPTION:
3783 * Given an array of 'num' double values, which have been
3784 * sorted, find the value corresponding to the value which is at
3785 * the 'perc'th percentile in the list array. Return this value.
3786 *
3787 * RETURN:
3788 * double value of the number at 'perc'th percentile in array
3789 *
3790 * </AUTO>
3791 */
3792
3793 static double
find_percentile(double * array,int num,double perc)3794 find_percentile
3795 (
3796 double *array, /* I: look in this SORTED array */
3797 int num, /* I: which has this many elements */
3798 double perc /* I: for entry at this percentile */
3799 )
3800 {
3801 int index;
3802
3803 shAssert(array != NULL);
3804 shAssert(num > 0);
3805 shAssert((perc > 0.0) && (perc <= 1.0));
3806
3807 index = (int) floor(num*perc + 0.5);
3808 if (index >= num) {
3809 index = num - 1;
3810 }
3811 return(array[index]);
3812 }
3813
3814
3815 /************************************************************************
3816 *
3817 *
3818 * ROUTINE: calc_trans_coords
3819 *
3820 * DESCRIPTION:
3821 * Given a single s_star structure, apply the
3822 * given TRANS structure to its coordinates.
3823 * Place the converted coordinates into the given output args
3824 * "newx" and "newy".
3825 *
3826 * We use the trans->order value to flag the type of transformation
3827 * to calculate.
3828 *
3829 *
3830 * RETURN:
3831 * SH_SUCCESS if all goes well
3832 * SH_GENERIC_ERROR if some problem occurs
3833 *
3834 * </AUTO>
3835 */
3836
3837 static int
calc_trans_coords(s_star * star,TRANS * trans,double * newx,double * newy)3838 calc_trans_coords
3839 (
3840 s_star *star, /* I: use this STAR's coords as input */
3841 TRANS *trans, /* I: contains coefficients of transformation */
3842 double *newx, /* O: contains output x coord */
3843 double *newy /* O: contains output y coord */
3844 )
3845 {
3846 double rsquared;
3847
3848 shAssert(star != NULL);
3849 shAssert(trans != NULL);
3850
3851 switch (trans->order) {
3852 case AT_TRANS_LINEAR:
3853 *newx = trans->a + trans->b*star->x + trans->c*star->y;
3854 *newy = trans->d + trans->e*star->x + trans->f*star->y;
3855 break;
3856
3857 case AT_TRANS_QUADRATIC:
3858 *newx = trans->a + trans->b*star->x + trans->c*star->y +
3859 trans->d*star->x*star->x + trans->e*star->x*star->y +
3860 trans->f*star->y*star->y;
3861 *newy = trans->g + trans->h*star->x + trans->i*star->y +
3862 trans->j*star->x*star->x + trans->k*star->x*star->y +
3863 trans->l*star->y*star->y;
3864 break;
3865
3866 case AT_TRANS_CUBIC:
3867 rsquared = star->x*star->x + star->y*star->y;
3868 *newx = trans->a + trans->b*star->x + trans->c*star->y +
3869 trans->d*star->x*star->x + trans->e*star->x*star->y +
3870 trans->f*star->y*star->y +
3871 trans->g*star->x*rsquared + trans->h*star->y*rsquared;
3872
3873 *newy = trans->i + trans->j*star->x + trans->k*star->y +
3874 trans->l*star->x*star->x + trans->m*star->x*star->y +
3875 trans->n*star->y*star->y +
3876 trans->o*star->x*rsquared + trans->p*star->y*rsquared;
3877 break;
3878
3879 default:
3880 shFatal("calc_trans_coords: given invalid trans->order %d \n",
3881 trans->order);
3882 break;
3883 }
3884
3885 return(SH_SUCCESS);
3886 }
3887
3888
3889 /************************************************************************
3890 *
3891 *
3892 * ROUTINE: apply_trans
3893 *
3894 * DESCRIPTION:
3895 * Given an array of 'num_stars' s_star structures, apply the
3896 * given TRANS structure to the coordinates of each one.
3897 *
3898 *
3899 * RETURN:
3900 * SH_SUCCESS if all goes well
3901 * SH_GENERIC_ERROR if some problem occurs
3902 *
3903 * </AUTO>
3904 */
3905
3906 static int
apply_trans(s_star * star_array,int num_stars,TRANS * trans)3907 apply_trans
3908 (
3909 s_star *star_array, /* I/O: array of structures to modify */
3910 int num_stars, /* I: number of stars in the array */
3911 TRANS *trans /* I: contains coefficients of transformation */
3912 )
3913 {
3914 int i;
3915 double newx, newy;
3916 s_star *sp;
3917
3918 if (num_stars == 0) {
3919 return(SH_SUCCESS);
3920 }
3921 shAssert(star_array != NULL);
3922 shAssert(trans != NULL);
3923
3924 for (i = 0; i < num_stars; i++) {
3925 sp = &(star_array[i]);
3926 if (calc_trans_coords(sp, trans, &newx, &newy) != SH_SUCCESS) {
3927 shError("apply_trans: calc_trans_coords fails");
3928 return(SH_GENERIC_ERROR);
3929 }
3930 sp->x = newx;
3931 sp->y = newy;
3932 }
3933
3934 return(SH_SUCCESS);
3935 }
3936
3937
3938
3939
3940
3941 /***************************************************************************
3942 *
3943 *
3944 * ROUTINE: double_sort_by_match_id
3945 *
3946 * DESCRIPTION:
3947 * sort all the elements of the first array of "s_star" in increasing
3948 * order by "match_id" value. Also, reorder the
3949 * elements of the _second_ array in exactly the same way, so that
3950 * the elements of both array which matched BEFORE the sorting
3951 * will match again _after_ the sorting.
3952 *
3953 * return:
3954 * SH_SUCCESS if all goes well
3955 * SH_GENERIC_ERROR if not
3956 *
3957 * </AUTO>
3958 */
3959
3960 static int
double_sort_by_match_id(s_star * star_array_A,int num_stars_A,s_star * star_array_B,int num_stars_B)3961 double_sort_by_match_id
3962 (
3963 s_star *star_array_A, /* I/O: array to be sorted */
3964 int num_stars_A, /* I: number of stars in array A */
3965 s_star *star_array_B, /* I/O: array to be re-ordered just as A */
3966 int num_stars_B /* I: number of stars in array B */
3967 )
3968 {
3969 int i;
3970 struct s_star *temp_array;
3971 struct s_star *sb, *stemp;
3972
3973 shAssert(num_stars_A == num_stars_B);
3974 if (num_stars_A == 0) {
3975 return(SH_SUCCESS);
3976 }
3977 shAssert(star_array_A != NULL);
3978 shAssert(star_array_B != NULL);
3979
3980 /*
3981 * first, let's set the "index" field of each element of each
3982 * star_array its position in the array.
3983 */
3984 for (i = 0; i < num_stars_A; i++) {
3985 star_array_A[i].index = i;
3986 star_array_B[i].index = i;
3987 }
3988
3989 /*
3990 * next, we create a temporary array of the same size as A and B.
3991 */
3992 temp_array = (s_star *) shMalloc(num_stars_A*sizeof(s_star));
3993
3994
3995 /*
3996 * Now, the two arrays A and B are currently arranged so that
3997 * star_array_A[i] matches star_array_B[i]. We want to sort
3998 * star_array_A, and re-arrange star_array_B so that the
3999 * corresponding elements still match up afterwards.
4000 *
4001 * - sort star_array_A
4002 * - loop i through sorted star_array_A
4003 * copy star_array_B element matching star_array_A[i]
4004 * into temp_array[i]
4005 * - loop i through star_array_B
4006 * copy temp_array[i] into star_array_B[i]
4007 *
4008 * - delete temp_array
4009 *
4010 * We end up with star_array_A sorted by "x", and star_array_B
4011 * re-arranged in exactly the same order.
4012 */
4013
4014 sort_star_by_match_id(star_array_A, num_stars_A);
4015 for (i = 0; i < num_stars_A; i++) {
4016 sb = &(star_array_B[star_array_A[i].index]);
4017 shAssert(sb != NULL);
4018 stemp = &(temp_array[i]);
4019 shAssert(stemp != NULL);
4020 copy_star(sb, stemp);
4021 }
4022
4023 /*
4024 * now copy the elements of the temp_array back into star_array_B
4025 */
4026 for (i = 0; i < num_stars_A; i++) {
4027 sb = &(star_array_B[i]);
4028 shAssert(sb != NULL);
4029 stemp = &(temp_array[i]);
4030 shAssert(stemp != NULL);
4031 copy_star(stemp, sb);
4032 }
4033
4034 /*
4035 * and we're done! Delete the temporary array
4036 */
4037 free_star_array(temp_array);
4038
4039 return(SH_SUCCESS);
4040 }
4041
4042
4043
4044
4045
4046
4047 /***************************************************************************
4048 *
4049 *
4050 * ROUTINE: match_arrays_slow
4051 *
4052 * DESCRIPTION:
4053 * given two arrays of s_stars [A and B], find all matching elements,
4054 * where a match is coincidence of centers to within "radius" pixels.
4055 *
4056 * Use a slow, but sure, algorithm (and an inefficient implementation,
4057 * I'm sure. As of 1/18/96, trying for correctness, not speed).
4058 *
4059 * We will match objects from A --> B. It is possible to have several
4060 * As that match to the same B:
4061 *
4062 * A1 -> B5 and A2 -> B5
4063 *
4064 * This function finds such multiple-match items and deletes all but
4065 * the closest of the matches.
4066 *
4067 * This array creates 4 new arrays of s_stars, and returns a pointer
4068 * to each array, as well as the number of stars in each array.
4069 *
4070 * place the elems of A that are matches into output array J
4071 * B that are matches into output array K
4072 * A that are not matches into output array L
4073 * B that are not matches into output array M
4074 *
4075 * return: SH_SUCCESS if all goes well
4076 * SH_GENERIC_ERROR if not
4077 *
4078 * </AUTO>
4079 */
4080
4081
4082 static int
match_arrays_slow(s_star * star_array_A,int num_stars_A,s_star * star_array_B,int num_stars_B,double radius,s_star ** star_array_J,int * num_stars_J,s_star ** star_array_K,int * num_stars_K,s_star ** star_array_L,int * num_stars_L,s_star ** star_array_M,int * num_stars_M)4083 match_arrays_slow
4084 (
4085 s_star *star_array_A, /* I: first array of s_stars to be matched */
4086 int num_stars_A, /* I: number of stars in A */
4087 s_star *star_array_B, /* I: second array of s_stars to be matched */
4088 int num_stars_B, /* I: number of stars in B */
4089 double radius, /* I: matching radius */
4090 s_star **star_array_J, /* O: all stars in A which match put in here */
4091 int *num_stars_J, /* O: number of stars in output array J */
4092 s_star **star_array_K, /* O: all stars in B which match put in here */
4093 int *num_stars_K, /* O: number of stars in output array K */
4094 s_star **star_array_L, /* O: all stars in A which don't match put here */
4095 int *num_stars_L, /* O: number of stars in output array L */
4096 s_star **star_array_M, /* O: all stars in B which don't match put here */
4097 int *num_stars_M /* O: number of stars in output array M */
4098 )
4099 {
4100 double Ax, Ay, Bx, By;
4101 double dist, limit;
4102 int matchPos;
4103 int posA, posB;
4104 int current_num_J, current_num_K;
4105 double deltax, deltay;
4106 double Axm, Axp, Aym, Ayp;
4107 s_star *sa, *sb;
4108
4109 #ifdef DEBUG
4110 printf("entering match_arrays_slow ");
4111 #endif
4112
4113 /*
4114 * first, we create each of the 4 output arrays. We start with
4115 * each as big as the input arrays, but we'll shrink them down
4116 * to their proper sizes before we return.
4117 */
4118 *star_array_J = (s_star *) shMalloc(num_stars_A*sizeof(s_star));
4119 *num_stars_J = num_stars_A;
4120 *star_array_K = (s_star *) shMalloc(num_stars_B*sizeof(s_star));
4121 *num_stars_K = num_stars_B;
4122 *star_array_L = (s_star *) shMalloc(num_stars_A*sizeof(s_star));
4123 *num_stars_L = num_stars_A;
4124 *star_array_M = (s_star *) shMalloc(num_stars_B*sizeof(s_star));
4125 *num_stars_M = num_stars_B;
4126
4127 /*
4128 * make some sanity checks
4129 */
4130 shAssert(num_stars_A >= 0);
4131 shAssert(num_stars_B >= 0);
4132 if ((num_stars_A == 0) || (num_stars_B == 0)) {
4133 return(SH_SUCCESS);
4134 }
4135 shAssert(star_array_A != NULL);
4136 shAssert(star_array_B != NULL);
4137
4138
4139 matchPos = 0; /* placate compiler */
4140
4141
4142 /*
4143 * First, we sort arrays A and B by their "x" coordinates,
4144 * to facilitate matching.
4145 */
4146 sort_star_by_x(star_array_A, num_stars_A);
4147 sort_star_by_x(star_array_B, num_stars_B);
4148
4149 /*
4150 * We copy array A into L, and array B into M.
4151 * We will remove all non-matching elements from these
4152 * output arrays later on in this function.
4153 */
4154
4155 copy_star_array(star_array_A, *star_array_L, num_stars_A);
4156 copy_star_array(star_array_B, *star_array_M, num_stars_B);
4157
4158
4159 /*
4160 * this is the largest distance that two stars can be from
4161 * each other and still be a match.
4162 */
4163 limit = radius*radius;
4164
4165
4166 /*
4167 * the first step is to go slowly through array A, checking against
4168 * every object in array B. If there's a match, we copy the matching
4169 * elements onto lists J and K, respectively. We do NOT check
4170 * yet to see if there are multiply-matched elements.
4171 *
4172 * This implementation could be speeded up a LOT by sorting the
4173 * two arrays in "x" and then making use of the information to check
4174 * only stars which are close to each other in "x". Do that
4175 * some time later.... MWR 1/18/96.
4176 */
4177 #ifdef DEBUG
4178 printf(" size of array A is %d, array B is %d\n", num_stars_A, num_stars_B);
4179 printf(" about to step through array A looking for matches\n");
4180 #endif
4181
4182 current_num_J = 0;
4183 current_num_K = 0;
4184
4185 for (posA = 0; posA < num_stars_A; posA++) {
4186
4187 shAssert((sa = &(star_array_A[posA])) != NULL);
4188 Ax = sa->x;
4189 Ay = sa->y;
4190
4191 Axm = Ax - radius;
4192 Axp = Ax + radius;
4193 Aym = Ay - radius;
4194 Ayp = Ay + radius;
4195
4196 for (posB = 0; posB < num_stars_B; posB++) {
4197
4198 shAssert((sb = &(star_array_B[posB])) != NULL);
4199 Bx = sb->x;
4200 By = sb->y;
4201
4202 /* check quickly to see if we can avoid a multiply */
4203 if ((Bx < Axm) || (Bx > Axp) || (By < Aym) || (By > Ayp)) {
4204 continue;
4205 }
4206
4207 /* okay, we actually have to calculate a distance here. */
4208 deltax = Ax - Bx;
4209 deltay = Ay - By;
4210 dist = deltax*deltax + deltay*deltay;
4211 if (dist < limit) {
4212
4213 /*
4214 * we have a match (at least, a possible match). So, copy
4215 * objA onto listJ and objB onto listK. But do NOT remove
4216 * these objects from listA and listB! We may end up
4217 * matching another objA to the same objB later on, and
4218 * we will continue trying to match this same objA to other
4219 * objBs.
4220 */
4221 add_element(sa, star_array_J, num_stars_J, ¤t_num_J);
4222 add_element(sb, star_array_K, num_stars_K, ¤t_num_K);
4223
4224 }
4225 }
4226 }
4227
4228 /*
4229 * at this point, let's re-set "*num_stars_J" to the proper number.
4230 * Recall that the "add_element" function may increase "*num_stars_J"
4231 * by factors of 2, while the variable "current_num_J" keeps track
4232 * of the actual number of stars in the array. It ought to be the
4233 * case that
4234 * num_stars_J <= *num_stars_J
4235 *
4236 * and likewise for K.
4237 */
4238 *num_stars_J = current_num_J;
4239 *num_stars_K = current_num_K;
4240
4241 #ifdef DEBUG
4242 printf(" done with stepping through array A \n");
4243 printf(" array J has %d, array K has %d \n", current_num_J, current_num_K);
4244 #endif
4245
4246 #ifdef DEBUG
4247 /* for debugging only */
4248 for (posA = 0; posA < *num_stars_J; posA++) {
4249 sa = &((*star_array_J)[posA]);
4250 sb = &((*star_array_K)[posA]);
4251 printf(" %4d J: %4d (%8.2f, %8.2f) K: %4d (%8.2f, %8.2f) \n",
4252 posA, sa->match_id, sa->x, sa->y, sb->match_id, sb->x, sb->y);
4253 }
4254 #endif
4255
4256
4257 /*
4258 * at this point, all _possible_ matches have been placed into
4259 * corresponding elements of arrays J and K. Now, we go through
4260 * array J to find elements which appear more than once. We'll
4261 * resolve them by throwing out all but the closest match.
4262 */
4263
4264 /*
4265 * first, sort array J by the "match_id" values. This allows us to find
4266 * repeated elements easily. Re-order array K in exactly the same
4267 * way, so matching elements still match.
4268 */
4269 #ifdef DEBUG
4270 printf(" sorting array J by match_id\n");
4271 #endif
4272 if (double_sort_by_match_id(*star_array_J, *num_stars_J,
4273 *star_array_K, *num_stars_K) != SH_SUCCESS) {
4274 shError("match_arrays_slow: can't sort array J");
4275 return(SH_GENERIC_ERROR);
4276 }
4277 #ifdef DEBUG
4278 for (posA = 0; posA < *num_stars_J; posA++) {
4279 sa = &((*star_array_J)[posA]);
4280 sb = &((*star_array_K)[posA]);
4281 printf(" %4d J: %4d (%8.2f, %8.2f) K: %4d (%8.2f, %8.2f) \n",
4282 posA, sa->match_id, sa->x, sa->y, sb->match_id, sb->x, sb->y);
4283 }
4284 #endif
4285
4286 /*
4287 * now remove repeated elements from array J, keeping the closest matches
4288 */
4289 #ifdef DEBUG
4290 printf(" before remove_repeated_elements, array J has %d\n", *num_stars_J);
4291 #endif
4292 if (remove_repeated_elements(*star_array_J, num_stars_J,
4293 *star_array_K, num_stars_K) != SH_SUCCESS) {
4294 shError("match_arrays_slow: remove_repeated_elements fails for array J");
4295 return(SH_GENERIC_ERROR);
4296 }
4297 #ifdef DEBUG
4298 printf(" after remove_repeated_elements, array J has %d\n", *num_stars_J);
4299 for (posA = 0; posA < *num_stars_J; posA++) {
4300 sa = &((*star_array_J)[posA]);
4301 sb = &((*star_array_K)[posA]);
4302 printf(" %4d J: %4d (%8.2f, %8.2f) K: %4d (%8.2f, %8.2f) \n",
4303 posA, sa->match_id, sa->x, sa->y, sb->match_id, sb->x, sb->y);
4304 }
4305 #endif
4306 shAssert(*num_stars_J == *num_stars_K);
4307
4308 /*
4309 * next, do the same for array K: sort it by "match_id"
4310 * (and re-arrange array J to match),
4311 * then find and remove any
4312 * repeated elements, keeping only the closest matches.
4313 */
4314 #ifdef DEBUG
4315 printf(" sorting array K by match_id\n");
4316 #endif
4317 if (double_sort_by_match_id(*star_array_K, *num_stars_K,
4318 *star_array_J, *num_stars_J) != SH_SUCCESS) {
4319 shError("match_arrays_slow: can't sort array K");
4320 return(SH_GENERIC_ERROR);
4321 }
4322 #ifdef DEBUG
4323 for (posA = 0; posA < *num_stars_J; posA++) {
4324 sa = &((*star_array_J)[posA]);
4325 sb = &((*star_array_K)[posA]);
4326 printf(" %4d J: %4d (%8.2f, %8.2f) K: %4d (%8.2f, %8.2f) \n",
4327 posA, sa->match_id, sa->x, sa->y, sb->match_id, sb->x, sb->y);
4328 }
4329 #endif
4330
4331
4332 #ifdef DEBUG
4333 printf(" before remove_repeated_elements, array K has %d\n", *num_stars_K);
4334 for (posA = 0; posA < *num_stars_J; posA++) {
4335 sa = &((*star_array_J)[posA]);
4336 sb = &((*star_array_K)[posA]);
4337 printf(" %4d J: %4d (%8.2f, %8.2f) K: %4d (%8.2f, %8.2f) \n",
4338 posA, sa->match_id, sa->x, sa->y, sb->match_id, sb->x, sb->y);
4339 }
4340 #endif
4341 if (remove_repeated_elements(*star_array_K, num_stars_K,
4342 *star_array_J, num_stars_J) != SH_SUCCESS) {
4343 shError("match_arrays_slow: remove_repeated_elements fails for array K");
4344 return(SH_GENERIC_ERROR);
4345 }
4346 #ifdef DEBUG
4347 printf(" after remove_repeated_elements, arrary K has %d\n", *num_stars_K);
4348 for (posA = 0; posA < *num_stars_J; posA++) {
4349 sa = &((*star_array_J)[posA]);
4350 sb = &((*star_array_K)[posA]);
4351 printf(" %4d J: %4d (%8.2f, %8.2f) K: %4d (%8.2f, %8.2f) \n",
4352 posA, sa->match_id, sa->x, sa->y, sb->match_id, sb->x, sb->y);
4353 }
4354 #endif
4355 shAssert(*num_stars_J == *num_stars_K);
4356
4357 /*
4358 * finally, we have unique set of closest-pair matching elements
4359 * in arrays J and K. Now we can remove any element from array L
4360 * which appears in array J, and remove any element from array M
4361 * which appears in array K. First, we'll sort arrays L and M
4362 * to make the comparisons easier.
4363 */
4364 #ifdef DEBUG
4365 printf(" sorting array L \n");
4366 #endif
4367 sort_star_by_match_id(*star_array_L, *num_stars_L);
4368 #ifdef DEBUG
4369 printf(" sorting array M \n");
4370 #endif
4371 sort_star_by_match_id(*star_array_M, *num_stars_M);
4372
4373 /*
4374 * Recall that array K is already sorted by "match_id", but that
4375 * we may have thrown J out of order when we forced it to follow
4376 * the sorting of K. So, first we'll sort J by "match_id",
4377 * (and re-order K match it), then we can remove repeated elements
4378 * from L easily.
4379 */
4380 #ifdef DEBUG
4381 printf(" sorting array J by match_id\n");
4382 #endif
4383 if (double_sort_by_match_id(*star_array_J, *num_stars_J,
4384 *star_array_K, *num_stars_K) != SH_SUCCESS) {
4385 shError("match_arrays_slow: can't sort array J");
4386 return(SH_GENERIC_ERROR);
4387 }
4388 #ifdef DEBUG
4389 printf(" after double_sort_by_match_id (J, K)\n");
4390 for (posA = 0; posA < *num_stars_J; posA++) {
4391 sa = &((*star_array_J)[posA]);
4392 sb = &((*star_array_K)[posA]);
4393 printf(" %4d J: %4d (%8.2f, %8.2f) K: %4d (%8.2f, %8.2f) \n",
4394 posA, sa->match_id, sa->x, sa->y, sb->match_id, sb->x, sb->y);
4395 }
4396 #endif
4397 /*
4398 * now remove elements from array L which appear in array J
4399 */
4400 #ifdef DEBUG
4401 printf(" before remove_same_elements, array L has %d\n", *num_stars_L);
4402 for (posA = 0; posA < *num_stars_L; posA++) {
4403 sa = &((*star_array_L)[posA]);
4404 sb = &((*star_array_M)[posA]);
4405 printf(" %4d L: %4d (%8.2f, %8.2f) M: %4d (%8.2f, %8.2f) \n",
4406 posA, sa->match_id, sa->x, sa->y, sb->match_id, sb->x, sb->y);
4407 }
4408 #endif
4409 remove_same_elements(*star_array_J, *num_stars_J,
4410 *star_array_L, num_stars_L);
4411 #ifdef DEBUG
4412 printf(" after remove_same_elements, array L has %d\n", *num_stars_L);
4413 for (posA = 0; posA < *num_stars_L; posA++) {
4414 sa = &((*star_array_L)[posA]);
4415 sb = &((*star_array_M)[posA]);
4416 printf(" %4d L: %4d (%8.2f, %8.2f) M: %4d (%8.2f, %8.2f) \n",
4417 posA, sa->match_id, sa->x, sa->y, sb->match_id, sb->x, sb->y);
4418 }
4419 #endif
4420
4421
4422 /*
4423 * Recall that we threw K out of order when we forced it to follow
4424 * the sorting of J. So, we'll sort K by "match_id",
4425 * (and re-order J match it), then we can remove repeated elements
4426 * from M easily.
4427 */
4428 #ifdef DEBUG
4429 printf(" sorting array K by match_id\n");
4430 #endif
4431 if (double_sort_by_match_id(*star_array_K, *num_stars_K,
4432 *star_array_J, *num_stars_J) != SH_SUCCESS) {
4433 shError("match_arrays_slow: can't sort array K");
4434 return(SH_GENERIC_ERROR);
4435 }
4436 #ifdef DEBUG
4437 printf(" after double_sort_by_match_id (K, J)\n");
4438 for (posA = 0; posA < *num_stars_J; posA++) {
4439 sa = &((*star_array_J)[posA]);
4440 sb = &((*star_array_K)[posA]);
4441 printf(" %4d J: %4d (%8.2f, %8.2f) K: %4d (%8.2f, %8.2f) \n",
4442 posA, sa->match_id, sa->x, sa->y, sb->match_id, sb->x, sb->y);
4443 }
4444 #endif
4445 /*
4446 * and remove elements from array M which appear in array K
4447 */
4448 #ifdef DEBUG
4449 printf(" before remove_same_elements, array M has %d\n", *num_stars_M);
4450 for (posA = 0; posA < *num_stars_L; posA++) {
4451 sa = &((*star_array_L)[posA]);
4452 sb = &((*star_array_M)[posA]);
4453 printf(" %4d L: %4d (%8.2f, %8.2f) M: %4d (%8.2f, %8.2f) \n",
4454 posA, sa->match_id, sa->x, sa->y, sb->match_id, sb->x, sb->y);
4455 }
4456 #endif
4457 remove_same_elements(*star_array_K, *num_stars_K,
4458 *star_array_M, num_stars_M);
4459 #ifdef DEBUG
4460 printf(" after remove_same_elements, array M has %d\n", *num_stars_M);
4461 for (posA = 0; posA < *num_stars_L; posA++) {
4462 sa = &((*star_array_L)[posA]);
4463 sb = &((*star_array_M)[posA]);
4464 printf(" %4d L: %4d (%8.2f, %8.2f) M: %4d (%8.2f, %8.2f) \n",
4465 posA, sa->match_id, sa->x, sa->y, sb->match_id, sb->x, sb->y);
4466 }
4467 #endif
4468
4469
4470 return(SH_SUCCESS);
4471 }
4472
4473
4474
4475 /**************************************************************************
4476 *
4477 *
4478 * ROUTINE: add_element
4479 *
4480 * DESCRIPTION:
4481 * We are given a pointer to s_star, an array of "total_num" s_stars,
4482 * and a count of the current number of s_stars set in the array.
4483 *
4484 * We want to copy the contents of the single star into
4485 * the "current_num"'th element of the array.
4486 *
4487 * If current_num < total_num, just perform copy,
4488 * increment current_num
4489 *
4490 * If current_num == total_num, we must allocate more space in array
4491 * allocate an array 2x as big as total_num
4492 * copy existing elements into new array
4493 * copy new element into new array
4494 * free old array
4495 * make old array pointer point to new array
4496 * increment current_num
4497 *
4498 * We could avoid all this by using linked lists, but I think
4499 * that we will only rarely have to increase the size of an array,
4500 * and never increase its size more than once. So this isn't so bad.
4501 *
4502 * RETURN: SH_SUCCESS if all goes well
4503 * SH_GENERIC_ERROR if not
4504 *
4505 */
4506
4507 static int
add_element(s_star * new_star,s_star ** star_array,int * total_num,int * current_num)4508 add_element
4509 (
4510 s_star *new_star, /* I: want to copy this into next slot in array */
4511 s_star **star_array, /* I/O: will copy into this array */
4512 /* if necessary, will allocate a new array, */
4513 /* copy entire contents into it, including */
4514 /* new_star, then free the old array */
4515 int *total_num, /* I/O: total number of stars allocated in */
4516 /* star_array. We may increase this if */
4517 /* we have to extend star_array */
4518 int *current_num /* I/O: current number of stars in star_array */
4519 /* which have been set. This number should */
4520 /* always increase by 1 if we succeed in */
4521 /* in adding the "new_star" */
4522 )
4523 {
4524 int num;
4525 s_star *new_array;
4526
4527 shAssert(new_star != NULL);
4528 shAssert((star_array != NULL) && (*star_array != NULL));
4529 shAssert((total_num != NULL) && (*total_num >= 0));
4530 shAssert((current_num != NULL) && (*current_num >= 0));
4531
4532
4533 /*
4534 * check for the easy case: if current_num < total_num, we can
4535 * just set star_array[current_num] and increment current_num.
4536 */
4537 if (*current_num < *total_num) {
4538 copy_star(new_star, &((*star_array)[*current_num]));
4539 (*current_num)++;
4540 }
4541 else if (*current_num == *total_num) {
4542
4543 /*
4544 * this is the tricky case, in which we have to allocate space
4545 * for a larger array, copy all existing elements, and then
4546 * copy over the new_star.
4547 */
4548 num = (*total_num)*2;
4549 new_array = (s_star *) shMalloc(num*sizeof(s_star));
4550 copy_star_array((*star_array), new_array, (*total_num));
4551 free_star_array(*star_array);
4552 *star_array = new_array;
4553 *total_num = num;
4554 copy_star(new_star, &((*star_array)[*current_num]));
4555 (*current_num)++;
4556 }
4557 else {
4558
4559 /*
4560 * this should never occur!
4561 */
4562 shAssert(0);
4563 }
4564
4565 return(SH_SUCCESS);
4566 }
4567
4568
4569
4570
4571
4572 /*********************************************************************
4573 *
4574 * ROUTINE: remove_repeated_elements
4575 *
4576 * DESCRIPTION:
4577 * step through the first array argument, star_array_1, checking for
4578 * successive elements which are the same. for each such pair, calculate the
4579 * distance between the matching elements of objects in arrays 1 and 2.
4580 * Throw the less-close pair out of the two array, modifying the number
4581 * of elements in each accordingly (and moving all other elements
4582 * up one place in the array).
4583 *
4584 * The two arrays must have the same number of elements,
4585 * and array 1 must already have been sorted by the "match_id" field.
4586 *
4587 * RETURN:
4588 * SH_SUCCESS if all goes well
4589 * SH_GENERIC_ERROR if something goes wrong
4590 *
4591 */
4592
4593 static int
remove_repeated_elements(s_star * star_array_1,int * num_stars_1,s_star * star_array_2,int * num_stars_2)4594 remove_repeated_elements
4595 (
4596 s_star *star_array_1, /* I/O: look in this array for repeats */
4597 int *num_stars_1, /* I/O: number of stars in array 1 */
4598 s_star *star_array_2, /* I/O: do to this array what we do to array 1 */
4599 int *num_stars_2 /* I/O: number of stars in array 2 */
4600 )
4601 {
4602 int pos1, pos2;
4603 double thisdist, lastdist;
4604 s_star *s1, *s2;
4605 s_star *last1, *last2;
4606
4607 shAssert(star_array_1 != NULL);
4608 shAssert(star_array_2 != NULL);
4609 shAssert(*num_stars_1 == *num_stars_2);
4610
4611 pos1 = 0;
4612 pos2 = 0;
4613
4614 last1 = NULL;
4615 last2 = NULL;
4616 while (pos1 < *num_stars_1) {
4617
4618 s1 = &(star_array_1[pos1]);
4619 s2 = &(star_array_2[pos2]);
4620 if ((s1 == NULL) || (s2 == NULL)) {
4621 shError("remove_repeated_elements: missing elem in array 1 or 2");
4622 return(SH_GENERIC_ERROR);
4623 }
4624
4625 if (last1 == NULL) {
4626 last1 = s1;
4627 last2 = s2;
4628 }
4629 else if (s1->match_id == last1->match_id) {
4630
4631 /* there is a repeated element. We must find the closer match */
4632 thisdist = (s1->x - s2->x)*(s1->x - s2->x) +
4633 (s1->y - s2->y)*(s1->y - s2->y);
4634 lastdist = (last1->x - last2->x)*(last1->x - last2->x) +
4635 (last1->y - last2->y)*(last1->y - last2->y);
4636
4637 if (thisdist < lastdist) {
4638
4639 /*
4640 * remove the "last" item from arrays 1 and 2.
4641 * We move the "current" items up one position in the arrays,
4642 * (into spaces [pos1 - 1] and [pos2 - 1]), and make
4643 * them the new "last" items.
4644 */
4645 remove_elem(star_array_1, pos1 - 1, num_stars_1);
4646 remove_elem(star_array_2, pos2 - 1, num_stars_2);
4647 last1 = &(star_array_1[pos1 - 1]);
4648 last2 = &(star_array_2[pos2 - 1]);
4649 }
4650 else {
4651
4652 /*
4653 * remove the current item from arrays 1 and 2.
4654 * We can leave the "last" items as they are, since
4655 * we haven't moved them.
4656 */
4657 remove_elem(star_array_1, pos1, num_stars_1);
4658 remove_elem(star_array_2, pos2, num_stars_2);
4659 }
4660 pos1--;
4661 pos2--;
4662 }
4663 else {
4664
4665 /* no repeated element. Prepare for next step forward */
4666 last1 = s1;
4667 last2 = s2;
4668 }
4669 pos1++;
4670 pos2++;
4671 }
4672 return(SH_SUCCESS);
4673 }
4674
4675
4676 /*********************************************************************
4677 *
4678 * ROUTINE: remove_elem
4679 *
4680 * DESCRIPTION:
4681 * Remove the i'th element from the given array.
4682 *
4683 * What we do (slow as it is) is
4684 *
4685 * 1. move all elements after i up by 1
4686 * 2. subtract 1 from the number of elements in the array
4687 *
4688 * There's probably a better way of doing this, but let's
4689 * go with it for now. 1/19/96 MWR
4690 *
4691 * RETURN:
4692 * nothing
4693 *
4694 */
4695
4696 static void
remove_elem(s_star * star_array,int num,int * num_stars)4697 remove_elem
4698 (
4699 s_star *star_array, /* I/O: we remove one element from this array */
4700 int num, /* I: remove _this_ element */
4701 int *num_stars /* I/O: on input: number of stars in array */
4702 /* on output: ditto, now smaller by one */
4703 )
4704 {
4705 int i;
4706 s_star *s1, *s2;
4707
4708 shAssert(star_array != NULL);
4709 shAssert(num < *num_stars);
4710 shAssert(num >= 0);
4711 shAssert(*num_stars > 0);
4712
4713 s1 = &(star_array[num]);
4714 s2 = &(star_array[num + 1]);
4715 for (i = num; i < ((*num_stars) - 1); i++, s1++, s2++) {
4716 copy_star(s2, s1);
4717 }
4718
4719 (*num_stars)--;
4720 }
4721
4722
4723 /*********************************************************************
4724 *
4725 * ROUTINE: remove_same_elements
4726 *
4727 * DESCRIPTION:
4728 * given two arrays of s_stars which have been sorted by their
4729 * "match_id" values, try to find s_stars which appear
4730 * in both arrays. Remove any such s_stars from the second array.
4731 *
4732 * RETURN:
4733 * nothing
4734 *
4735 */
4736
4737 static void
remove_same_elements(s_star * star_array_1,int num_stars_1,s_star * star_array_2,int * num_stars_2)4738 remove_same_elements
4739 (
4740 s_star *star_array_1, /* I: look for elems which match those in array 2 */
4741 int num_stars_1, /* I: number of elems in array 1 */
4742 s_star *star_array_2, /* I/O: remove elems which match those in array 1 */
4743 int *num_stars_2 /* I/O: number of elems in array 2 */
4744 /* will probably be smaller on output */
4745 )
4746 {
4747 int pos1, pos2, pos2_top;
4748 s_star *s1, *s2;
4749
4750 shAssert(star_array_1 != NULL);
4751 shAssert(star_array_2 != NULL);
4752 shAssert(num_stars_2 != NULL);
4753
4754 pos1 = 0;
4755 pos2_top = 0;
4756
4757 while (pos1 < num_stars_1) {
4758
4759 s1 = &(star_array_1[pos1]);
4760 shAssert(s1 != NULL);
4761
4762 for (pos2 = pos2_top; pos2 < *num_stars_2; pos2++) {
4763 s2 = &(star_array_2[pos2]);
4764 shAssert(s2 != NULL);
4765
4766 if (s1->match_id == s2->match_id) {
4767 remove_elem(star_array_2, pos2, num_stars_2);
4768 if (--pos2_top < 0) {
4769 pos2_top = 0;
4770 }
4771 }
4772 else {
4773 if (s2->match_id < s1->match_id) {
4774 pos2_top = pos2 + 1;
4775 }
4776 }
4777 }
4778 pos1++;
4779 }
4780 }
4781
4782
4783 /***********************************************************************
4784 * ROUTINE: list_to_array
4785 *
4786 * DESCRIPTION:
4787 * Create an array of s_star structures, identical to the given linked
4788 * list. Just make a copy of each structure.
4789 *
4790 * Sure, this is inefficient, but I'm using legacy code ...
4791 *
4792 * Return a pointer to the complete, filled array.
4793 *
4794 */
4795
4796 static s_star *
list_to_array(int num_stars,struct s_star * list)4797 list_to_array
4798 (
4799 int num_stars, /* I: number of stars in the list */
4800 struct s_star *list /* I: the linked list */
4801 )
4802 {
4803 int i;
4804 struct s_star *array = NULL;
4805 struct s_star *ptr;
4806 struct s_star *star;
4807
4808 /*
4809 * okay, now we can walk down the CHAIN and create a new s_star
4810 * for each item on the CHAIN.
4811 */
4812 array = (s_star *) shMalloc(num_stars*sizeof(s_star));
4813 shAssert(array != NULL);
4814 for (i = 0, ptr = list; i < num_stars; i++, ptr = ptr->next) {
4815 shAssert(ptr != NULL);
4816 star = &(array[i]);
4817 shAssert(star != NULL);
4818 set_star(star, ptr->x, ptr->y, ptr->mag);
4819 star->match_id = i;
4820 }
4821
4822 return(array);
4823 }
4824
4825
4826 /***********************************************************************
4827 * ROUTINE: write_array
4828 *
4829 * DESCRIPTION:
4830 * Given an array of s_star structures, write them to an ASCII text
4831 * file, with the following format:
4832 *
4833 * ID xvalue yvalue magvalue
4834 *
4835 * The 'ID' value is one assigned internally by these routines --
4836 * it doesn't correspond to any input ID value.
4837 *
4838 * RETURNS:
4839 * nothing
4840 */
4841
4842 static void
write_array(int num_stars,struct s_star * star_array,char * filename)4843 write_array
4844 (
4845 int num_stars, /* I: number of stars in the array */
4846 struct s_star *star_array, /* I: the array of stars */
4847 char *filename /* I: write into this file */
4848 )
4849 {
4850 int i;
4851 FILE *fp;
4852
4853 if ((fp = fopen(filename, "w")) == NULL) {
4854 shFatal("write_array: can't open file %s", filename);
4855 }
4856
4857 for (i = 0; i < num_stars; i++) {
4858 fprintf(fp, "%6d %13.7f %13.7f %8.4f\n", star_array[i].id,
4859 star_array[i].x, star_array[i].y, star_array[i].mag);
4860 }
4861
4862 fclose(fp);
4863 }
4864
4865
4866 /****************************************************************************
4867 * ROUTINE: write_small_arrays
4868 *
4869 * Given an array of 'nstar' stars, and an array of s_triangle structures,
4870 * write them into a file which is a mixture of ASCII and binary data.
4871 * We use a format which is ASCII at first, but turns to binary at the end:
4872 *
4873 * ASCII section:
4874 * line 1 nstar
4875 * line 2 info on star 1
4876 * line 3 info on star 2
4877 * ...
4878 * line nstar+1 info on star nstar
4879 * line nstar+2 ntriangle
4880 *
4881 * BINARY section:
4882 * binary dump of triangle 1
4883 * binary dump of triangle 2
4884 * ...
4885 * binary dump of triangle ntriangle
4886 *
4887 *
4888 * RETURN:
4889 * SH_SUCCESS if all goes well
4890 * SH_GENERIC_ERORR if not
4891 */
4892
4893 static int
write_small_arrays(double ra,double dec,int num_stars,struct s_star * star_array,int nbright,int num_triangles,struct s_triangle * t_array,char * outfile)4894 write_small_arrays
4895 (
4896 double ra, /* I: Right Ascension of field, degrees */
4897 double dec, /* I: Declination of field, degrees */
4898 int num_stars, /* I: number of stars in linked list */
4899 struct s_star *star_array, /* I: array of s_star structures */
4900 int nbright, /* I: write only this many stars, at most */
4901 int num_triangles, /* I: number of triangles in array */
4902 struct s_triangle *t_array, /* I: array of triangle structures */
4903 char *outfile /* I: name of output file */
4904 )
4905 {
4906 int i, num;
4907 struct s_star *star;
4908 struct s_triangle *tri;
4909 FILE *fp;
4910
4911 if (num_stars > 0) {
4912 shAssert(star_array != NULL);
4913 }
4914 if (num_triangles > 0) {
4915 shAssert(t_array != NULL);
4916 }
4917
4918 num = (num_stars < nbright ? num_stars : nbright);
4919
4920 if ((fp = fopen(outfile, "w")) == NULL) {
4921 shError("write_small_arrays: can't open file %s for writing", outfile);
4922 return(SH_GENERIC_ERROR);
4923 }
4924
4925 /* write out the RA and Dec of field center, in decimal degrees */
4926 fprintf(fp, "%f %f\n", ra, dec);
4927
4928 /* write out the stars */
4929 fprintf(fp, "%d\n", num);
4930 for (i = 0; i < num; i++) {
4931 star = &(star_array[i]);
4932 shAssert(star != NULL);
4933
4934 fprintf(fp, "%6d %6d %12.5f %12.5f %12.5f\n",
4935 i, star->id, star->x, star->y, star->mag);
4936 }
4937
4938 /* write out the triangles */
4939 fprintf(fp, "%d\n", num_triangles);
4940 for (i = 0; i < num_triangles; i++) {
4941 tri = &(t_array[i]);
4942 shAssert(tri != NULL);
4943
4944 #ifdef MATCH_ASCII
4945 /* this old version creates an ASCII file, which takes forever to read */
4946 fprintf(fp, "%6d %6d %13.5f %6.4f %6.4f %4d %4d %4d\n",
4947 i, tri->id, tri->a_length,
4948 tri->ba, tri->ca,
4949 tri->a_index,
4950 tri->b_index,
4951 tri->c_index);
4952 #else
4953 /* but this binary version can be read in a flash */
4954 tri->index = i;
4955 tri->match_id = -1;
4956 tri->next = NULL;
4957 fwrite((void *) tri, sizeof(struct s_triangle), 1, fp);
4958 #endif
4959
4960 }
4961 fclose(fp);
4962
4963 return(SH_SUCCESS);
4964 }
4965
4966
4967
4968 /***********************************************************************
4969 * FUNCTION: reset_array_ids
4970 *
4971 * Modify the 'id' field values in the given array
4972 * so that they will match the 'id' values in the corresponding
4973 * stars of the given list.
4974 *
4975 * RETURNS
4976 * nothing
4977 */
4978
4979 static void
reset_array_ids(struct s_star * star_list,int num_stars,struct s_star * star_array)4980 reset_array_ids
4981 (
4982 struct s_star *star_list, /* I: a list of stars */
4983 int num_stars, /* I: number of stars in list and array */
4984 struct s_star *star_array /* I/O: reset 'id' fields in this array */
4985 )
4986 {
4987 int i;
4988 struct s_star *star_in_list, *star_in_array;
4989
4990 star_in_list = star_list;
4991 for (i = 0; i < num_stars; i++) {
4992
4993 star_in_array = &(star_array[i]);
4994 shAssert(star_in_list != NULL);
4995 shAssert(star_in_array != NULL);
4996
4997 star_in_array->id = star_in_list->id;
4998
4999 star_in_list = star_in_list->next;
5000 }
5001 }
5002
5003
5004
5005
5006
5007 /************************************************************************
5008 *
5009 *
5010 * ROUTINE: calc_trans_linear
5011 *
5012 * DESCRIPTION:
5013 * Given a set of "nbright" matched pairs of stars, which we can
5014 * extract from the "winner_index" and "star_array" arrays,
5015 * figure out a TRANS structure which takes coordinates of
5016 * objects in set A and transforms then into coords for set B.
5017 *
5018 * In this case, a TRANS contains 6 coefficients in equations like this:
5019 *
5020 * x' = A + B*x + C*y
5021 * y' = D + E*x + F*y
5022 *
5023 * where (x,y) are coords in set A and (x',y') are corresponding
5024 * coords in set B.
5025 *
5026 * Internally, I'm going to solve for the very similar equations
5027 *
5028 * x' = Ax + By + C
5029 * y' = Dx + Ey + F
5030 *
5031 * and then just re-arrange the coefficients at the very end. OK?
5032 *
5033 *
5034 * What we do is to treat each of the two equations above
5035 * separately. We can write down 3 equations relating quantities
5036 * in the two sets of points (there are more than 3 such equations,
5037 * but we don't seek an exhaustive list). For example,
5038 *
5039 * a. x' = Ax + By + C
5040 * b. x'x = Ax^2 + Bxy + Cx (mult both sides by x)
5041 * c. x'y = Axy + By^2 + Cy (mult both sides by y)
5042 *
5043 * Now, since we have "nbright" matched pairs, we can take each of
5044 * the above 3 equations and form the sums on both sides, over
5045 * all "nbright" points. So, if S(x) represents the sum of the quantity
5046 * "x" over all nbright points, and if we let N=nbright, then
5047 *
5048 * a. S(x') = AS(x) + BS(y) + CN
5049 * b. S(x'x) = AS(x^2) + BS(xy) + CS(x)
5050 * c. S(x'y) = AS(xy) + BS(y^2) + CS(y)
5051 *
5052 * At this point, we have a set of three equations, and 3 unknowns: A, B, C.
5053 * We can write this set of equations as a matrix equation
5054 *
5055 * b = M * v
5056 *
5057 * where we KNOW the quantities
5058 *
5059 * vector b = ( S(x'), S(x'x), S(x'y) )
5060 *
5061 * matrix M = ( S(x) S(y) 1 )
5062 * ( S(x^2) S(xy) S(x) )
5063 * ( S(xy) S(y^2) S(y) )
5064 *
5065 *
5066 * and we want to FIND the unknown
5067 *
5068 * vector v = ( A, B, C )
5069 *
5070 * So, how to solve this matrix equation? We use a Gaussian-elimination
5071 * method (see notes in 'gauss_matrix' function). We solve
5072 * for A, B, C (and equivalently for D, E, F), then fill in the fields
5073 * of the given TRANS structure argument.
5074 *
5075 * It's possible that the matrix will be singular, and we can't find
5076 * a solution. In that case, we print an error message and don't touch
5077 * the TRANS' fields.
5078 *
5079 * [should explain how we make an iterative solution here,
5080 * but will put in comments later. MWR ]
5081 *
5082 * RETURN:
5083 * SH_SUCCESS if all goes well
5084 * SH_GENERIC_ERROR if we can't find a solution
5085 *
5086 * </AUTO>
5087 */
5088
5089 static int
calc_trans_linear(int nbright,s_star * star_array_A,int num_stars_A,s_star * star_array_B,int num_stars_B,int * winner_votes,int * winner_index_A,int * winner_index_B,TRANS * trans)5090 calc_trans_linear
5091 (
5092 int nbright, /* I: max number of stars we use in calculating */
5093 /* the transformation; we may cut down to */
5094 /* a more well-behaved subset. */
5095 s_star *star_array_A, /* I: first array of s_star structure we match */
5096 /* the output TRANS takes their coords */
5097 /* into those of array B */
5098 int num_stars_A, /* I: total number of stars in star_array_A */
5099 s_star *star_array_B, /* I: second array of s_star structure we match */
5100 int num_stars_B, /* I: total number of stars in star_array_B */
5101 int *winner_votes, /* I: number of votes gotten by the top 'nbright' */
5102 /* matched pairs of stars */
5103 int *winner_index_A, /* I: index into "star_array_A" of top */
5104 /* vote-getters */
5105 int *winner_index_B, /* I: index into "star_array_B" of top */
5106 /* vote-getters */
5107 TRANS *trans /* O: place solved coefficients into this */
5108 /* existing structure's fields */
5109 )
5110 {
5111 int i;
5112 double **matrix;
5113 double vector[3];
5114 double solved_a, solved_b, solved_c, solved_d, solved_e, solved_f;
5115 s_star *s1, *s2;
5116 /* */
5117 double sum, sumx1, sumy1, sumx2, sumy2;
5118 double sumx1sq, sumy1sq;
5119 double sumx1y1, sumx1x2, sumx1y2;
5120 double sumy1x2, sumy1y2;
5121
5122
5123 shAssert(nbright >= AT_MATCH_REQUIRE_LINEAR);
5124 shAssert(trans->order == AT_TRANS_LINEAR);
5125
5126
5127 /*
5128 * allocate a matrix we'll need for this function
5129 */
5130 matrix = alloc_matrix(3);
5131
5132
5133 /*
5134 * first, we consider the coefficients A, B, C in the trans.
5135 * we form the sums that make up the elements of matrix M
5136 */
5137 sum = 0.0;
5138 sumx1 = 0.0;
5139 sumy1 = 0.0;
5140 sumx2 = 0.0;
5141 sumy2 = 0.0;
5142 sumx1sq = 0.0;
5143 sumy1sq = 0.0;
5144 sumx1x2 = 0.0;
5145 sumx1y1 = 0.0;
5146 sumx1y2 = 0.0;
5147 sumy1x2 = 0.0;
5148 sumy1y2 = 0.0;
5149
5150 for (i = 0; i < nbright; i++) {
5151
5152 /* sanity checks */
5153 shAssert(winner_index_A[i] < num_stars_A);
5154 s1 = &(star_array_A[winner_index_A[i]]);
5155 shAssert(winner_index_B[i] < num_stars_B);
5156 s2 = &(star_array_B[winner_index_B[i]]);
5157
5158 /* elements of the matrix */
5159 sum += 1.0;
5160 sumx1 += s1->x;
5161 sumx2 += s2->x;
5162 sumy1 += s1->y;
5163 sumy2 += s2->y;
5164 sumx1sq += s1->x*s1->x;
5165 sumy1sq += s1->y*s1->y;
5166 sumx1x2 += s1->x*s2->x;
5167 sumx1y1 += s1->x*s1->y;
5168 sumx1y2 += s1->x*s2->y;
5169 sumy1x2 += s1->y*s2->x;
5170 sumy1y2 += s1->y*s2->y;
5171
5172 }
5173
5174
5175 /*
5176 * now turn these sums into a matrix and a vector
5177 */
5178 matrix[0][0] = sumx1sq;
5179 matrix[0][1] = sumx1y1;
5180 matrix[0][2] = sumx1;
5181 matrix[1][0] = sumx1y1;
5182 matrix[1][1] = sumy1sq;
5183 matrix[1][2] = sumy1;
5184 matrix[2][0] = sumx1;
5185 matrix[2][1] = sumy1;
5186 matrix[2][2] = sum;
5187
5188 vector[0] = sumx1x2;
5189 vector[1] = sumy1x2;
5190 vector[2] = sumx2;
5191
5192 #ifdef DEBUG
5193 printf("before calling solution routines for ABC, here's matrix\n");
5194 print_matrix(matrix, 3);
5195 #endif
5196
5197 /*
5198 * and now call the Gaussian-elimination routines to solve the matrix.
5199 * The solution for TRANS coefficients A, B, C will be placed
5200 * into the elements on "vector" after "gauss_matrix" finishes.
5201 */
5202 if (gauss_matrix(matrix, 3, vector) != SH_SUCCESS) {
5203 shError("calc_trans_linear: can't solve for coeffs A,B,C ");
5204 return(SH_GENERIC_ERROR);
5205 }
5206
5207 #ifdef DEBUG
5208 printf("after calling solution routines, here's matrix\n");
5209 print_matrix(matrix, 3);
5210 #endif
5211
5212 solved_a = vector[0];
5213 solved_b = vector[1];
5214 solved_c = vector[2];
5215
5216
5217 /*
5218 * Okay, now we solve for TRANS coefficients D, E, F, using the
5219 * set of equations that relates y' to (x,y)
5220 *
5221 * a. y' = Dx + Ey + F
5222 * b. y'x = Dx^2 + Exy + Fx (mult both sides by x)
5223 * c. y'y = Dxy + Ey^2 + Fy (mult both sides by y)
5224 *
5225 */
5226 matrix[0][0] = sumx1sq;
5227 matrix[0][1] = sumx1y1;
5228 matrix[0][2] = sumx1;
5229 matrix[1][0] = sumx1y1;
5230 matrix[1][1] = sumy1sq;
5231 matrix[1][2] = sumy1;
5232 matrix[2][0] = sumx1;
5233 matrix[2][1] = sumy1;
5234 matrix[2][2] = sum;
5235
5236 vector[0] = sumx1y2;
5237 vector[1] = sumy1y2;
5238 vector[2] = sumy2;
5239
5240 #ifdef DEBUG
5241 printf("before calling solution routines for DEF, here's matrix\n");
5242 print_matrix(matrix, 3);
5243 #endif
5244
5245 /*
5246 * and now call the Gaussian-elimination routines to solve the matrix.
5247 * The solution for TRANS coefficients D, E, F will be placed
5248 * into the elements on "vector" after "gauss_matrix" finishes.
5249 */
5250 if (gauss_matrix(matrix, 3, vector) != SH_SUCCESS) {
5251 shError("calc_trans_linear: can't solve for coeffs D,E,F ");
5252 return(SH_GENERIC_ERROR);
5253 }
5254
5255 #ifdef DEBUG
5256 printf("after calling solution routines, here's matrix\n");
5257 print_matrix(matrix, 3);
5258 #endif
5259
5260 solved_d = vector[0];
5261 solved_e = vector[1];
5262 solved_f = vector[2];
5263
5264
5265 /*
5266 * assign the coefficients we've just calculated to the output
5267 * TRANS structure. Recall that we've solved equations
5268 *
5269 * x' = Ax + By + C
5270 * y' = Dx + Ey + F
5271 *
5272 * but that the TRANS structure assigns its coefficients assuming
5273 *
5274 * x' = A + Bx + Cy
5275 * y' = D + Ex + Fy
5276 *
5277 * so, here, we have to re-arrange the coefficients a bit.
5278 */
5279 trans->a = solved_c;
5280 trans->b = solved_a;
5281 trans->c = solved_b;
5282 trans->d = solved_f;
5283 trans->e = solved_d;
5284 trans->f = solved_e;
5285
5286 /*
5287 * free up memory we allocated for this function
5288 */
5289 free_matrix(matrix, 3);
5290
5291 return(SH_SUCCESS);
5292 }
5293
5294
5295 /************************************************************************
5296 *
5297 *
5298 * ROUTINE: calc_trans_quadratic
5299 *
5300 * DESCRIPTION:
5301 * Given a set of "nbright" matched pairs of stars, which we can
5302 * extract from the "winner_index" and "star_array" arrays,
5303 * figure out a TRANS structure which takes coordinates of
5304 * objects in set A and transforms then into coords for set B.
5305 * In this case, a TRANS contains the twelve coefficients in the equations
5306 *
5307 * x' = A + Bx + Cy + Dxx + Exy + Fyy
5308 * y' = G + Hx + Iy + Jxx + Kxy + Lyy
5309 *
5310 * where (x,y) are coords in set A and (x',y') are corresponding
5311 * coords in set B.
5312 *
5313 *
5314 * What we do is to treat each of the two equations above
5315 * separately. We can write down 6 equations relating quantities
5316 * in the two sets of points (there are more than 6 such equations,
5317 * but we don't seek an exhaustive list). For example,
5318 *
5319 * a. x' = A + Bx + Cy + Dxx + Exy + Fyy
5320 * b. x'x = Ax + Bxx + Cxy + Dxxx + Exxy + Fxyy
5321 * c. x'y = Ay + Bxy + Cyy + Dxxy + Exyy + Fyyy
5322 * d. x'xx = Axx + Bxxx + Cxxy + Dxxxx + Exxxy + Fxxyy
5323 * e. x'xy = Axy + Bxxy + Cxyy + Dxxxy + Exxyy + Fxyyy
5324 * f. x'yy = Ayy + Bxyy + Cyyy + Dxxyy + Exyyy + Fyyyy
5325 *
5326 * Now, since we have "nbright" matched pairs, we can take each of
5327 * the above 6 equations and form the sums on both sides, over
5328 * all "nbright" points. So, if S(x) represents the sum of the quantity
5329 * "x" over all nbright points, and if we let N=nbright, then
5330 *
5331 * a. S(x') = AN + BS(x) + CS(y) + DS(xx) + ES(xy) + FS(yy)
5332 * b. S(x'x) = AS(x) + BS(xx) + CS(xy) + DS(xxx) + ES(xxy) + FS(xyy)
5333 * c. S(x'y) = AS(y) + BS(xy) + CS(yy) + DS(xxy) + ES(xyy) + FS(yyy)
5334 * d. S(x'xx) = AS(xx) + BS(xxx) + CS(xxy) + DS(xxxx) + ES(xxxy) + FS(xxyy)
5335 * e. S(x'xy) = AS(xy) + BS(xxy) + CS(xyy) + DS(xxxy) + ES(xxyy) + FS(xyyy)
5336 * f. S(x'yy) = AS(yy) + BS(xyy) + CS(yyy) + DS(xxyy) + ES(xyyy) + FS(yyyy)
5337 *
5338 * At this point, we have a set of 6 equations, and 6 unknowns:
5339 * A, B, C, D, E, F
5340 *
5341 * We can write this set of equations as a matrix equation
5342 *
5343 * b = M * v
5344 *
5345 * where we KNOW the quantities
5346 *
5347 * vector b = ( S(x'), S(x'x), S(x'y), S(x'xx), S(x'xy), S(x'yy) )
5348 *
5349 * matrix M = [ N S(x) S(y) S(xx) S(xy) S(yy) ]
5350 * [ S(x) S(xx) S(xy) S(xxx) S(xxy) S(xyy) ]
5351 * [ S(y) S(xy) S(yy) S(xxy) S(xyy) S(yyy) ]
5352 * [ S(xx) S(xxx) S(xxy) S(xxxx) S(xxxy) S(xxyy) ]
5353 * [ S(xy) S(xxy) S(xyy) S(xxxy) S(xxyy) S(xyyy) ]
5354 * [ S(yy) S(xyy) S(yyy) S(xxyy) S(xyyy) S(yyyy) ]
5355 *
5356 * and we want to FIND the unknown
5357 *
5358 * vector v = ( A, B, C, D, E, F )
5359 *
5360 * So, how to solve this matrix equation? We use a Gaussian-elimination
5361 * method (see notes in 'gauss_matrix' function). We solve
5362 * for A, B, C, D, E, F (and equivalently for G, H, I, J, K, L),
5363 * then fill in the fields
5364 * of the given TRANS structure argument.
5365 *
5366 * It's possible that the matrix will be singular, and we can't find
5367 * a solution. In that case, we print an error message and don't touch
5368 * the TRANS' fields.
5369 *
5370 * [should explain how we make an iterative solution here,
5371 * but will put in comments later. MWR ]
5372 *
5373 * RETURN:
5374 * SH_SUCCESS if all goes well
5375 * SH_GENERIC_ERROR if we can't find a solution
5376 *
5377 * </AUTO>
5378 */
5379
5380 static int
calc_trans_quadratic(int nbright,s_star * star_array_A,int num_stars_A,s_star * star_array_B,int num_stars_B,int * winner_votes,int * winner_index_A,int * winner_index_B,TRANS * trans)5381 calc_trans_quadratic
5382 (
5383 int nbright, /* I: max number of stars we use in calculating */
5384 /* the transformation; we may cut down to */
5385 /* a more well-behaved subset. */
5386 s_star *star_array_A, /* I: first array of s_star structure we match */
5387 /* the output TRANS takes their coords */
5388 /* into those of array B */
5389 int num_stars_A, /* I: total number of stars in star_array_A */
5390 s_star *star_array_B, /* I: second array of s_star structure we match */
5391 int num_stars_B, /* I: total number of stars in star_array_B */
5392 int *winner_votes, /* I: number of votes gotten by the top 'nbright' */
5393 /* matched pairs of stars */
5394 int *winner_index_A, /* I: index into "star_array_A" of top */
5395 /* vote-getters */
5396 int *winner_index_B, /* I: index into "star_array_B" of top */
5397 /* vote-getters */
5398 TRANS *trans /* O: place solved coefficients into this */
5399 /* existing structure's fields */
5400 )
5401 {
5402 int i;
5403 double **matrix;
5404 double vector[6];
5405 double solved_a, solved_b, solved_c, solved_d, solved_e, solved_f;
5406 double solved_g, solved_h, solved_i, solved_j, solved_k, solved_l;
5407 s_star *s1, *s2;
5408
5409 /*
5410 * in variable names below, a '1' refers to coordinate of star s1
5411 * (which appear on both sides of the matrix equation)
5412 * and a '2' refers to coordinate of star s2
5413 * (which appears only on left hand side of matrix equation) o
5414 */
5415 double sumx2, sumx2x1, sumx2y1, sumx2x1sq, sumx2x1y1, sumx2y1sq;
5416 double sumy2, sumy2x1, sumy2y1, sumy2x1sq, sumy2x1y1, sumy2y1sq;
5417
5418 double sum, sumx1, sumy1, sumx1sq, sumx1y1, sumy1sq;
5419 double sumx1cu, sumx1sqy1, sumx1y1sq;
5420 double sumy1cu;
5421 double sumx1qu, sumx1cuy1, sumx1sqy1sq;
5422 double sumx1y1cu;
5423 double sumy1qu;
5424
5425
5426 shAssert(nbright >= AT_MATCH_REQUIRE_QUADRATIC);
5427 shAssert(trans->order == AT_TRANS_QUADRATIC);
5428
5429
5430 /*
5431 * allocate a matrix we'll need for this function
5432 */
5433 matrix = alloc_matrix(6);
5434
5435
5436 /*
5437 * first, we consider the coefficients A, B, C, D, E, F in the trans.
5438 * we form the sums that make up the elements of matrix M
5439 */
5440
5441 sum = 0.0;
5442 sumx1 = 0.0;
5443 sumy1 = 0.0;
5444 sumx1sq = 0.0;
5445 sumx1y1 = 0.0;
5446 sumy1sq = 0.0;
5447 sumx1cu = 0.0;
5448 sumx1sqy1 = 0.0;
5449 sumx1y1sq = 0.0;
5450 sumy1cu = 0.0;
5451 sumx1qu = 0.0;
5452 sumx1cuy1 = 0.0;
5453 sumx1sqy1sq = 0.0;
5454 sumx1y1cu = 0.0;
5455 sumy1qu = 0.0;
5456
5457 sumx2 = 0.0; sumx2x1 = 0.0; sumx2y1 = 0.0;
5458 sumx2x1sq = 0.0; sumx2x1y1 = 0.0; sumx2y1sq = 0.0;
5459 sumy2 = 0.0; sumy2x1 = 0.0; sumy2y1 = 0.0;
5460 sumy2x1sq = 0.0; sumy2x1y1 = 0.0; sumy2y1sq = 0.0;
5461
5462
5463 for (i = 0; i < nbright; i++) {
5464
5465 /* sanity checks */
5466 shAssert(winner_index_A[i] < num_stars_A);
5467 s1 = &(star_array_A[winner_index_A[i]]);
5468 shAssert(winner_index_B[i] < num_stars_B);
5469 s2 = &(star_array_B[winner_index_B[i]]);
5470
5471 /* elements of the vectors */
5472 sumx2 += s2->x;
5473 sumx2x1 += s2->x*s1->x;
5474 sumx2y1 += s2->x*s1->y;
5475 sumx2x1sq += s2->x*s1->x*s1->x;
5476 sumx2x1y1 += s2->x*s1->x*s1->y;
5477 sumx2y1sq += s2->x*s1->y*s1->y;
5478
5479 sumy2 += s2->y;
5480 sumy2x1 += s2->y*s1->x;
5481 sumy2y1 += s2->y*s1->y;
5482 sumy2x1sq += s2->y*s1->x*s1->x;
5483 sumy2x1y1 += s2->y*s1->x*s1->y;
5484 sumy2y1sq += s2->y*s1->y*s1->y;
5485
5486
5487 /* elements of the matrix */
5488 sum += 1.0;
5489 sumx1 += s1->x;
5490 sumy1 += s1->y;
5491
5492 sumx1sq += s1->x*s1->x;
5493 sumx1y1 += s1->x*s1->y;
5494 sumy1sq += s1->y*s1->y;
5495
5496 sumx1cu += s1->x*s1->x*s1->x;
5497 sumx1sqy1 += s1->x*s1->x*s1->y;
5498 sumx1y1sq += s1->x*s1->y*s1->y;
5499 sumy1cu += s1->y*s1->y*s1->y;
5500
5501 sumx1qu += s1->x*s1->x*s1->x*s1->x;
5502 sumx1cuy1 += s1->x*s1->x*s1->x*s1->y;
5503 sumx1sqy1sq += s1->x*s1->x*s1->y*s1->y;
5504 sumx1y1cu += s1->x*s1->y*s1->y*s1->y;
5505 sumy1qu += s1->y*s1->y*s1->y*s1->y;
5506
5507 }
5508
5509
5510 /*
5511 * now turn these sums into a matrix and a vector
5512 */
5513 matrix[0][0] = sum;
5514 matrix[0][1] = sumx1;
5515 matrix[0][2] = sumy1;
5516 matrix[0][3] = sumx1sq;
5517 matrix[0][4] = sumx1y1;
5518 matrix[0][5] = sumy1sq;
5519
5520 matrix[1][0] = sumx1;
5521 matrix[1][1] = sumx1sq;
5522 matrix[1][2] = sumx1y1;
5523 matrix[1][3] = sumx1cu;
5524 matrix[1][4] = sumx1sqy1;
5525 matrix[1][5] = sumx1y1sq;
5526
5527 matrix[2][0] = sumy1;
5528 matrix[2][1] = sumx1y1;
5529 matrix[2][2] = sumy1sq;
5530 matrix[2][3] = sumx1sqy1;
5531 matrix[2][4] = sumx1y1sq;
5532 matrix[2][5] = sumy1cu;
5533
5534 matrix[3][0] = sumx1sq;
5535 matrix[3][1] = sumx1cu;
5536 matrix[3][2] = sumx1sqy1;
5537 matrix[3][3] = sumx1qu;
5538 matrix[3][4] = sumx1cuy1;
5539 matrix[3][5] = sumx1sqy1sq;
5540
5541 matrix[4][0] = sumx1y1;
5542 matrix[4][1] = sumx1sqy1;
5543 matrix[4][2] = sumx1y1sq;
5544 matrix[4][3] = sumx1cuy1;
5545 matrix[4][4] = sumx1sqy1sq;
5546 matrix[4][5] = sumx1y1cu;
5547
5548 matrix[5][0] = sumy1sq;
5549 matrix[5][1] = sumx1y1sq;
5550 matrix[5][2] = sumy1cu;
5551 matrix[5][3] = sumx1sqy1sq;
5552 matrix[5][4] = sumx1y1cu;
5553 matrix[5][5] = sumy1qu;
5554
5555 vector[0] = sumx2;
5556 vector[1] = sumx2x1;
5557 vector[2] = sumx2y1;
5558 vector[3] = sumx2x1sq;
5559 vector[4] = sumx2x1y1;
5560 vector[5] = sumx2y1sq;
5561
5562 #ifdef DEBUG
5563 printf("before calling solution routines for ABCDEF, here's matrix\n");
5564 print_matrix(matrix, 6);
5565 #endif
5566
5567 /*
5568 * and now call the Gaussian-elimination routines to solve the matrix.
5569 * The solution for TRANS coefficients A, B, C, D, E, F will be placed
5570 * into the elements on "vector" after "gauss_matrix" finishes.
5571 */
5572 if (gauss_matrix(matrix, 6, vector) != SH_SUCCESS) {
5573 shError("calc_trans_quadratic: can't solve for coeffs A,B,C,D,E,F ");
5574 return(SH_GENERIC_ERROR);
5575 }
5576
5577 #ifdef DEBUG
5578 printf("after calling solution routines, here's matrix\n");
5579 print_matrix(matrix, 6);
5580 #endif
5581
5582 solved_a = vector[0];
5583 solved_b = vector[1];
5584 solved_c = vector[2];
5585 solved_d = vector[3];
5586 solved_e = vector[4];
5587 solved_f = vector[5];
5588
5589
5590 /*
5591 * Okay, now we solve for TRANS coefficients G, H, I, J, K, L, using the
5592 * set of equations that relates y' to (x,y)
5593 *
5594 * y' = G + Hx + Iy + Jxx + Kxy + Lyy
5595 * y'x = Gx + Hxx + Ixy + Jxxx + Kxxy + Lxyy
5596 * y'y = Gy + Hxy + Iyy + Jxxy + Kxyy + Lyyy
5597 * y'xx = Gxx + Hxxx + Ixxy + Jxxxx + Kxxxy + Lxxyy
5598 * y'xy = Gxy + Hxxy + Ixyy + Jxxxy + Kxxyy + Lxyyy
5599 * y'yy = Gyy + Hxyy + Iyyy + Jxxyy + Kxyyy + Lyyyy
5600 *
5601 */
5602 matrix[0][0] = sum;
5603 matrix[0][1] = sumx1;
5604 matrix[0][2] = sumy1;
5605 matrix[0][3] = sumx1sq;
5606 matrix[0][4] = sumx1y1;
5607 matrix[0][5] = sumy1sq;
5608
5609 matrix[1][0] = sumx1;
5610 matrix[1][1] = sumx1sq;
5611 matrix[1][2] = sumx1y1;
5612 matrix[1][3] = sumx1cu;
5613 matrix[1][4] = sumx1sqy1;
5614 matrix[1][5] = sumx1y1sq;
5615
5616 matrix[2][0] = sumy1;
5617 matrix[2][1] = sumx1y1;
5618 matrix[2][2] = sumy1sq;
5619 matrix[2][3] = sumx1sqy1;
5620 matrix[2][4] = sumx1y1sq;
5621 matrix[2][5] = sumy1cu;
5622
5623 matrix[3][0] = sumx1sq;
5624 matrix[3][1] = sumx1cu;
5625 matrix[3][2] = sumx1sqy1;
5626 matrix[3][3] = sumx1qu;
5627 matrix[3][4] = sumx1cuy1;
5628 matrix[3][5] = sumx1sqy1sq;
5629
5630 matrix[4][0] = sumx1y1;
5631 matrix[4][1] = sumx1sqy1;
5632 matrix[4][2] = sumx1y1sq;
5633 matrix[4][3] = sumx1cuy1;
5634 matrix[4][4] = sumx1sqy1sq;
5635 matrix[4][5] = sumx1y1cu;
5636
5637 matrix[5][0] = sumy1sq;
5638 matrix[5][1] = sumx1y1sq;
5639 matrix[5][2] = sumy1cu;
5640 matrix[5][3] = sumx1sqy1sq;
5641 matrix[5][4] = sumx1y1cu;
5642 matrix[5][5] = sumy1qu;
5643
5644 vector[0] = sumy2;
5645 vector[1] = sumy2x1;
5646 vector[2] = sumy2y1;
5647 vector[3] = sumy2x1sq;
5648 vector[4] = sumy2x1y1;
5649 vector[5] = sumy2y1sq;
5650
5651 #ifdef DEBUG
5652 printf("before calling solution routines for GHIJKL, here's matrix\n");
5653 print_matrix(matrix, 6);
5654 #endif
5655
5656 /*
5657 * and now call the Gaussian-elimination routines to solve the matrix.
5658 * The solution for TRANS coefficients G, H, I, J, K, L will be placed
5659 * into the elements on "vector" after "gauss_matrix" finishes.
5660 */
5661 if (gauss_matrix(matrix, 6, vector) != SH_SUCCESS) {
5662 shError("calc_trans_quadratic: can't solve for coeffs G,H,I,J,K,L ");
5663 return(SH_GENERIC_ERROR);
5664 }
5665
5666 #ifdef DEBUG
5667 printf("after calling solution routines, here's matrix\n");
5668 print_matrix(matrix, 6);
5669 #endif
5670
5671 solved_g = vector[0];
5672 solved_h = vector[1];
5673 solved_i = vector[2];
5674 solved_j = vector[3];
5675 solved_k = vector[4];
5676 solved_l = vector[5];
5677
5678
5679 /*
5680 * assign the coefficients we've just calculated to the output
5681 * TRANS structure.
5682 */
5683 trans->a = solved_a;
5684 trans->b = solved_b;
5685 trans->c = solved_c;
5686 trans->d = solved_d;
5687 trans->e = solved_e;
5688 trans->f = solved_f;
5689 trans->g = solved_g;
5690 trans->h = solved_h;
5691 trans->i = solved_i;
5692 trans->j = solved_j;
5693 trans->k = solved_k;
5694 trans->l = solved_l;
5695
5696 /*
5697 * free up memory we allocated for this function
5698 */
5699 free_matrix(matrix, 6);
5700
5701 return(SH_SUCCESS);
5702 }
5703
5704
5705 /************************************************************************
5706 *
5707 *
5708 * ROUTINE: calc_trans_cubic
5709 *
5710 * DESCRIPTION:
5711 * Given a set of "nbright" matched pairs of stars, which we can
5712 * extract from the "winner_index" and "star_array" arrays,
5713 * figure out a TRANS structure which takes coordinates of
5714 * objects in set A and transforms then into coords for set B.
5715 * In this case, a TRANS contains the sixteen coefficients in the equations
5716 *
5717 * x' = A + Bx + Cy + Dxx + Exy + Fyy + Gx(xx+yy) + Hy(xx+yy)
5718 * y' = I + Jx + Ky + Lxx + Mxy + Nyy + Ox(xx+yy) + Py(xx+yy)
5719 *
5720 * where (x,y) are coords in set A and (x',y') are corresponding
5721 * coords in set B.
5722 *
5723 *
5724 * What we do is to treat each of the two equations above
5725 * separately. We can write down 8 equations relating quantities
5726 * in the two sets of points (there are more than 8 such equations,
5727 * but we don't seek an exhaustive list). For example,
5728 *
5729 * x' = A + Bx + Cy + Dxx + Exy + Fyy + GxR + HyR
5730 * x'x = Ax + Bxx + Cxy + Dxxx + Exxy + Fxyy + GxxR + HxyR
5731 * x'y = Ay + Bxy + Cyy + Dxxy + Exyy + Fyyy + GxyR + HyyR
5732 * x'xx = Axx + Bxxx + Cxxy + Dxxxx + Exxxy + Fxxyy + GxxxR + HxxyR
5733 * x'xy = Axy + Bxxy + Cxyy + Dxxxy + Exxyy + Fxyyy + GxxyR + HxyyR
5734 * x'yy = Ayy + Bxyy + Cyyy + Dxxyy + Exyyy + Fyyyy + GxyyR + HyyyR
5735 * x'xR = AxR + BxxR + CxyR + DxxxR + ExxyR + FxyyR + GxxRR + HxyRR
5736 * x'yR = AyR + BxyR + CyyR + DxxyR + ExyyR + FyyyR + GxyRR + HyyRR
5737 *
5738 * (where we have used 'R' as an abbreviation for (xx + yy))
5739 *
5740 * Now, since we have "nbright" matched pairs, we can take each of
5741 * the above 8 equations and form the sums on both sides, over
5742 * all "nbright" points. So, if S(x) represents the sum of the quantity
5743 * "x" over all nbright points, and if we let N=nbright, then
5744 *
5745 * S(x') = AN + BS(x) + CS(y) + DS(xx) + ES(xy) + FS(yy)
5746 * + GS(xR) + HS(yR)
5747 * S(x'x) = AS(x) + BS(xx) + CS(xy) + DS(xxx) + ES(xxy) + FS(xyy)
5748 * + GS(xxR) + HS(xyR)
5749 * S(x'y) = AS(y) + BS(xy) + CS(yy) + DS(xxy) + ES(xyy) + FS(yyy)
5750 * + GS(xyR) + HS(yyR)
5751 * S(x'xx) = AS(xx) + BS(xxx) + CS(xxy) + DS(xxxx) + ES(xxxy) + FS(xxyy)
5752 * + GS(xxxR) + HS(xxyR)
5753 * S(x'xy) = AS(xy) + BS(xxy) + CS(xyy) + DS(xxxy) + ES(xxyy) + FS(xyyy)
5754 * + GS(xxyR) + HS(xyyR)
5755 * S(x'yy) = AS(yy) + BS(xyy) + CS(yyy) + DS(xxyy) + ES(xyyy) + FS(yyyy)
5756 * + GS(xyyR) + HS(yyyR)
5757 * S(x'xR) = AS(xR) + BS(xxR) + CS(xyR) + DS(xxxR) + ES(xxyR) + FS(xyyR)
5758 * + GS(xxRR) + HS(xyRR)
5759 * S(x'yR) = AS(yR) + BS(xyR) + CS(yyR) + DS(xxyR) + ES(xyyR) + FS(yyyR)
5760 * + GS(xyRR) + HS(yyRR)
5761 *
5762 * At this point, we have a set of 8 equations, and 8 unknowns:
5763 * A, B, C, D, E, F, G, H
5764 *
5765 * We can write this set of equations as a matrix equation
5766 *
5767 * b = M * v
5768 *
5769 * where we KNOW the quantities
5770 *
5771 * b = ( S(x'), S(x'x), S(x'y), S(x'xx), S(x'xy), S(x'yy), S(x'xR), S(x'rR) )
5772 *
5773 * matr M = [ N S(x) S(y) S(xx) S(xy) S(yy) S(xR) S(yR) ]
5774 * [ S(x) S(xx) S(xy) S(xxx) S(xxy) S(xyy) S(xxR) S(xyR) ]
5775 * [ S(y) S(xy) S(yy) S(xxy) S(xyy) S(yyy) S(xyR) S(yyR) ]
5776 * [ S(xx) S(xxx) S(xxy) S(xxxx) S(xxxy) S(xxyy) S(xxxR) S(xxyR) ]
5777 * [ S(xy) S(xxy) S(xyy) S(xxxy) S(xxyy) S(xyyy) S(xxyR) S(xyyR) ]
5778 * [ S(yy) S(xyy) S(yyy) S(xxyy) S(xyyy) S(yyyy) S(xyyR) S(yyyR) ]
5779 * [ S(xR) S(xxR) S(xyR) S(xxxR) S(xxyR) S(xyyR) S(xxRR) S(xyRR) ]
5780 * [ S(yR) S(xyR) S(yyR) S(xxyR) S(xyyR) S(yyyR) S(xyRR) S(yyRR) ]
5781 *
5782 * and we want to FIND the unknown
5783 *
5784 * vector v = ( A, B, C, D, E, F, G, H )
5785 *
5786 * So, how to solve this matrix equation? We use a Gaussian-elimination
5787 * method (see notes in 'gauss_matrix' function). We solve
5788 * for A, B, C, D, E, F, G, H (and equivalently for I, J, K, L, M, N, O, P),
5789 * then fill in the fields
5790 * of the given TRANS structure argument.
5791 *
5792 * It's possible that the matrix will be singular, and we can't find
5793 * a solution. In that case, we print an error message and don't touch
5794 * the TRANS' fields.
5795 *
5796 * [should explain how we make an iterative solution here,
5797 * but will put in comments later. MWR ]
5798 *
5799 * RETURN:
5800 * SH_SUCCESS if all goes well
5801 * SH_GENERIC_ERROR if we can't find a solution
5802 *
5803 * </AUTO>
5804 */
5805
5806 static int
calc_trans_cubic(int nbright,s_star * star_array_A,int num_stars_A,s_star * star_array_B,int num_stars_B,int * winner_votes,int * winner_index_A,int * winner_index_B,TRANS * trans)5807 calc_trans_cubic
5808 (
5809 int nbright, /* I: max number of stars we use in calculating */
5810 /* the transformation; we may cut down to */
5811 /* a more well-behaved subset. */
5812 s_star *star_array_A, /* I: first array of s_star structure we match */
5813 /* the output TRANS takes their coords */
5814 /* into those of array B */
5815 int num_stars_A, /* I: total number of stars in star_array_A */
5816 s_star *star_array_B, /* I: second array of s_star structure we match */
5817 int num_stars_B, /* I: total number of stars in star_array_B */
5818 int *winner_votes, /* I: number of votes gotten by the top 'nbright' */
5819 /* matched pairs of stars */
5820 int *winner_index_A, /* I: index into "star_array_A" of top */
5821 /* vote-getters */
5822 int *winner_index_B, /* I: index into "star_array_B" of top */
5823 /* vote-getters */
5824 TRANS *trans /* O: place solved coefficients into this */
5825 /* existing structure's fields */
5826 )
5827 {
5828 int i;
5829 double **matrix;
5830 double vector[8];
5831 double solved_a, solved_b, solved_c, solved_d, solved_e, solved_f;
5832 double solved_g, solved_h;
5833 double solved_i, solved_j, solved_k, solved_l, solved_m, solved_n;
5834 double solved_o, solved_p;
5835 s_star *s1, *s2;
5836
5837 /*
5838 * the variable 'R' will hold the value (x1*x1 + y1*y1);
5839 * in other words, the square of the distance of (x1, y1)
5840 * from the origin.
5841 */
5842 double R;
5843 /*
5844 * in variable names below, a '1' refers to coordinate of star s1
5845 * (which appear on both sides of the matrix equation)
5846 * and a '2' refers to coordinate of star s2
5847 * (which appears only on left hand side of matrix equation) o
5848 */
5849 double sumx2, sumx2x1, sumx2y1, sumx2x1sq, sumx2x1y1, sumx2y1sq;
5850 double sumx2x1R, sumx2y1R;
5851 double sumy2, sumy2x1, sumy2y1, sumy2x1sq, sumy2x1y1, sumy2y1sq;
5852 double sumy2x1R, sumy2y1R;
5853
5854 double sum, sumx1, sumy1, sumx1sq, sumx1y1, sumy1sq;
5855 double sumx1cu, sumx1sqy1, sumx1y1sq;
5856 double sumy1cu;
5857 double sumx1R, sumy1R, sumx1sqR, sumx1y1R, sumy1sqR;
5858 double sumx1cuR, sumx1sqy1R, sumx1y1sqR, sumy1cuR;
5859 double sumx1qu, sumx1cuy1, sumx1sqy1sq;
5860 double sumx1y1cu;
5861 double sumy1qu;
5862 double sumx1sqRsq, sumx1y1Rsq, sumy1sqRsq;
5863
5864
5865 shAssert(nbright >= AT_MATCH_REQUIRE_CUBIC);
5866 shAssert(trans->order == AT_TRANS_CUBIC);
5867
5868
5869 /*
5870 * allocate a matrix we'll need for this function
5871 */
5872 matrix = alloc_matrix(8);
5873
5874
5875 /*
5876 * first, we consider the coefficients A, B, C, D, E, F, G, H in the trans.
5877 * we form the sums that make up the elements of matrix M
5878 */
5879
5880 sum = 0.0;
5881 sumx1 = 0.0;
5882 sumy1 = 0.0;
5883 sumx1sq = 0.0;
5884 sumx1y1 = 0.0;
5885 sumy1sq = 0.0;
5886 sumx1cu = 0.0;
5887 sumx1sqy1 = 0.0;
5888 sumx1y1sq = 0.0;
5889 sumy1cu = 0.0;
5890 sumx1qu = 0.0;
5891 sumx1cuy1 = 0.0;
5892 sumx1sqy1sq = 0.0;
5893 sumx1y1cu = 0.0;
5894 sumy1qu = 0.0;
5895 sumx1R = 0.0; sumy1R = 0.0; sumx1sqR = 0.0; sumx1y1R = 0.0; sumy1sqR = 0.0;
5896 sumx1cuR = 0.0; sumx1sqy1R = 0.0; sumx1y1sqR = 0.0; sumy1cuR = 0.0;
5897 sumx1sqRsq = 0.0; sumx1y1Rsq = 0.0; sumy1sqRsq = 0.0;
5898
5899 sumx2 = 0.0; sumx2x1 = 0.0; sumx2y1 = 0.0;
5900 sumx2x1sq = 0.0; sumx2x1y1 = 0.0; sumx2y1sq = 0.0;
5901 sumx2x1R = 0.0; sumx2y1R = 0.0;
5902 sumy2 = 0.0; sumy2x1 = 0.0; sumy2y1 = 0.0;
5903 sumy2x1sq = 0.0; sumy2x1y1 = 0.0; sumy2y1sq = 0.0;
5904 sumy2x1R = 0.0; sumy2y1R = 0.0;
5905
5906
5907 for (i = 0; i < nbright; i++) {
5908
5909 /* sanity checks */
5910 shAssert(winner_index_A[i] < num_stars_A);
5911 s1 = &(star_array_A[winner_index_A[i]]);
5912 shAssert(winner_index_B[i] < num_stars_B);
5913 s2 = &(star_array_B[winner_index_B[i]]);
5914
5915 /* elements of the vectors */
5916 R = (s1->x*s1->x + s1->y*s1->y);
5917
5918 sumx2 += s2->x;
5919 sumx2x1 += s2->x*s1->x;
5920 sumx2y1 += s2->x*s1->y;
5921 sumx2x1sq += s2->x*s1->x*s1->x;
5922 sumx2x1y1 += s2->x*s1->x*s1->y;
5923 sumx2y1sq += s2->x*s1->y*s1->y;
5924 sumx2x1R += s2->x*s1->x*R;
5925 sumx2y1R += s2->x*s1->y*R;
5926
5927 sumy2 += s2->y;
5928 sumy2x1 += s2->y*s1->x;
5929 sumy2y1 += s2->y*s1->y;
5930 sumy2x1sq += s2->y*s1->x*s1->x;
5931 sumy2x1y1 += s2->y*s1->x*s1->y;
5932 sumy2y1sq += s2->y*s1->y*s1->y;
5933 sumy2x1R += s2->y*s1->x*R;
5934 sumy2y1R += s2->y*s1->y*R;
5935
5936
5937 /* elements of the matrix */
5938 sum += 1.0;
5939 sumx1 += s1->x;
5940 sumy1 += s1->y;
5941
5942 sumx1sq += s1->x*s1->x;
5943 sumx1y1 += s1->x*s1->y;
5944 sumy1sq += s1->y*s1->y;
5945
5946 sumx1cu += s1->x*s1->x*s1->x;
5947 sumx1sqy1 += s1->x*s1->x*s1->y;
5948 sumx1y1sq += s1->x*s1->y*s1->y;
5949 sumy1cu += s1->y*s1->y*s1->y;
5950
5951 sumx1qu += s1->x*s1->x*s1->x*s1->x;
5952 sumx1cuy1 += s1->x*s1->x*s1->x*s1->y;
5953 sumx1sqy1sq += s1->x*s1->x*s1->y*s1->y;
5954 sumx1y1cu += s1->x*s1->y*s1->y*s1->y;
5955 sumy1qu += s1->y*s1->y*s1->y*s1->y;
5956
5957 sumx1R += s1->x*R;
5958 sumy1R += s1->y*R;
5959 sumx1sqR += s1->x*s1->x*R;
5960 sumx1y1R += s1->x*s1->y*R;
5961 sumy1sqR += s1->y*s1->y*R;
5962
5963 sumx1cuR += s1->x*s1->x*s1->x*R;
5964 sumx1sqy1R += s1->x*s1->x*s1->y*R;
5965 sumx1y1sqR += s1->x*s1->y*s1->y*R;
5966 sumy1cuR += s1->y*s1->y*s1->y*R;
5967
5968 sumx1sqRsq += s1->x*s1->x*R*R;
5969 sumx1y1Rsq += s1->x*s1->y*R*R;
5970 sumy1sqRsq += s1->y*s1->y*R*R;
5971
5972 }
5973
5974
5975 /*
5976 * now turn these sums into a matrix and a vector
5977 */
5978 matrix[0][0] = sum;
5979 matrix[0][1] = sumx1;
5980 matrix[0][2] = sumy1;
5981 matrix[0][3] = sumx1sq;
5982 matrix[0][4] = sumx1y1;
5983 matrix[0][5] = sumy1sq;
5984 matrix[0][6] = sumx1R;
5985 matrix[0][7] = sumy1R;
5986
5987 matrix[1][0] = sumx1;
5988 matrix[1][1] = sumx1sq;
5989 matrix[1][2] = sumx1y1;
5990 matrix[1][3] = sumx1cu;
5991 matrix[1][4] = sumx1sqy1;
5992 matrix[1][5] = sumx1y1sq;
5993 matrix[1][6] = sumx1sqR;
5994 matrix[1][7] = sumx1y1R;
5995
5996 matrix[2][0] = sumy1;
5997 matrix[2][1] = sumx1y1;
5998 matrix[2][2] = sumy1sq;
5999 matrix[2][3] = sumx1sqy1;
6000 matrix[2][4] = sumx1y1sq;
6001 matrix[2][5] = sumy1cu;
6002 matrix[2][6] = sumx1y1R;
6003 matrix[2][7] = sumy1sqR;
6004
6005 matrix[3][0] = sumx1sq;
6006 matrix[3][1] = sumx1cu;
6007 matrix[3][2] = sumx1sqy1;
6008 matrix[3][3] = sumx1qu;
6009 matrix[3][4] = sumx1cuy1;
6010 matrix[3][5] = sumx1sqy1sq;
6011 matrix[3][6] = sumx1cuR;
6012 matrix[3][7] = sumx1sqy1R;
6013
6014 matrix[4][0] = sumx1y1;
6015 matrix[4][1] = sumx1sqy1;
6016 matrix[4][2] = sumx1y1sq;
6017 matrix[4][3] = sumx1cuy1;
6018 matrix[4][4] = sumx1sqy1sq;
6019 matrix[4][5] = sumx1y1cu;
6020 matrix[4][6] = sumx1sqy1R;
6021 matrix[4][7] = sumx1y1sqR;
6022
6023 matrix[5][0] = sumy1sq;
6024 matrix[5][1] = sumx1y1sq;
6025 matrix[5][2] = sumy1cu;
6026 matrix[5][3] = sumx1sqy1sq;
6027 matrix[5][4] = sumx1y1cu;
6028 matrix[5][5] = sumy1qu;
6029 matrix[5][6] = sumx1y1sqR;
6030 matrix[5][7] = sumy1cuR;
6031
6032 matrix[6][0] = sumx1R;
6033 matrix[6][1] = sumx1sqR;
6034 matrix[6][2] = sumx1y1R;
6035 matrix[6][3] = sumx1cuR;
6036 matrix[6][4] = sumx1sqy1R;
6037 matrix[6][5] = sumx1y1sqR;
6038 matrix[6][6] = sumx1sqRsq;
6039 matrix[6][7] = sumx1y1Rsq;
6040
6041 matrix[7][0] = sumy1R;
6042 matrix[7][1] = sumx1y1R;
6043 matrix[7][2] = sumy1sqR;
6044 matrix[7][3] = sumx1sqy1R;
6045 matrix[7][4] = sumx1y1sqR;
6046 matrix[7][5] = sumy1cuR;
6047 matrix[7][6] = sumx1y1Rsq;
6048 matrix[7][7] = sumy1sqRsq;
6049
6050
6051 vector[0] = sumx2;
6052 vector[1] = sumx2x1;
6053 vector[2] = sumx2y1;
6054 vector[3] = sumx2x1sq;
6055 vector[4] = sumx2x1y1;
6056 vector[5] = sumx2y1sq;
6057 vector[6] = sumx2x1R;
6058 vector[7] = sumx2y1R;
6059
6060 #ifdef DEBUG
6061 printf("before calling solution routines for ABCDEFGH, here's matrix\n");
6062 print_matrix(matrix, 8);
6063 #endif
6064
6065 /*
6066 * and now call the Gaussian-elimination routines to solve the matrix.
6067 * The solution for TRANS coefficients A, B, C, D, E, F will be placed
6068 * into the elements on "vector" after "gauss_matrix" finishes.
6069 */
6070 if (gauss_matrix(matrix, 8, vector) != SH_SUCCESS) {
6071 shError("calc_trans_cubic: can't solve for coeffs A,B,C,D,E,F,G,H ");
6072 return(SH_GENERIC_ERROR);
6073 }
6074
6075 #ifdef DEBUG
6076 printf("after calling solution routines, here's matrix\n");
6077 print_matrix(matrix, 8);
6078 #endif
6079
6080 solved_a = vector[0];
6081 solved_b = vector[1];
6082 solved_c = vector[2];
6083 solved_d = vector[3];
6084 solved_e = vector[4];
6085 solved_f = vector[5];
6086 solved_g = vector[6];
6087 solved_h = vector[7];
6088
6089
6090 /*
6091 * Okay, now we solve for TRANS coefficients I, J, K, L, M, N, O, P
6092 * using the * set of equations that relates y' to (x,y)
6093 *
6094 * y' = I + Jx + Ky + Lxx + Mxy + Nyy + OxR + PyR
6095 * y'x = Ix + Jxx + Kxy + Lxxx + Mxxy + Nxyy + OxxR + PxyR
6096 * y'y = Iy + Jxy + Kyy + Lxxy + Mxyy + Nyyy + OxyR + PyyR
6097 * y'xx = Ixx + Jxxx + Kxxy + Lxxxx + Mxxxy + Nxxyy + OxxxR + PxxyR
6098 * y'xy = Ixy + Jxxy + Kxyy + Lxxxy + Mxxyy + Nxyyy + OxxyR + PxyyR
6099 * y'yy = Iyy + Jxyy + Kyyy + Lxxyy + Mxyyy + Nyyyy + OxyyR + PyyyR
6100 * y'xR = IxR + JxxR + KxyR + LxxxR + MxxyR + NxyyR + OxxRR + PxyRR
6101 * y'yR = IyR + JxyR + KyyR + LxxyR + MxyyR + NyyyR + OxyRR + PyyRR
6102 *
6103 */
6104 matrix[0][0] = sum;
6105 matrix[0][1] = sumx1;
6106 matrix[0][2] = sumy1;
6107 matrix[0][3] = sumx1sq;
6108 matrix[0][4] = sumx1y1;
6109 matrix[0][5] = sumy1sq;
6110 matrix[0][6] = sumx1R;
6111 matrix[0][7] = sumy1R;
6112
6113 matrix[1][0] = sumx1;
6114 matrix[1][1] = sumx1sq;
6115 matrix[1][2] = sumx1y1;
6116 matrix[1][3] = sumx1cu;
6117 matrix[1][4] = sumx1sqy1;
6118 matrix[1][5] = sumx1y1sq;
6119 matrix[1][6] = sumx1sqR;
6120 matrix[1][7] = sumx1y1R;
6121
6122 matrix[2][0] = sumy1;
6123 matrix[2][1] = sumx1y1;
6124 matrix[2][2] = sumy1sq;
6125 matrix[2][3] = sumx1sqy1;
6126 matrix[2][4] = sumx1y1sq;
6127 matrix[2][5] = sumy1cu;
6128 matrix[2][6] = sumx1y1R;
6129 matrix[2][7] = sumy1sqR;
6130
6131 matrix[3][0] = sumx1sq;
6132 matrix[3][1] = sumx1cu;
6133 matrix[3][2] = sumx1sqy1;
6134 matrix[3][3] = sumx1qu;
6135 matrix[3][4] = sumx1cuy1;
6136 matrix[3][5] = sumx1sqy1sq;
6137 matrix[3][6] = sumx1cuR;
6138 matrix[3][7] = sumx1sqy1R;
6139
6140 matrix[4][0] = sumx1y1;
6141 matrix[4][1] = sumx1sqy1;
6142 matrix[4][2] = sumx1y1sq;
6143 matrix[4][3] = sumx1cuy1;
6144 matrix[4][4] = sumx1sqy1sq;
6145 matrix[4][5] = sumx1y1cu;
6146 matrix[4][6] = sumx1sqy1R;
6147 matrix[4][7] = sumx1y1sqR;
6148
6149 matrix[5][0] = sumy1sq;
6150 matrix[5][1] = sumx1y1sq;
6151 matrix[5][2] = sumy1cu;
6152 matrix[5][3] = sumx1sqy1sq;
6153 matrix[5][4] = sumx1y1cu;
6154 matrix[5][5] = sumy1qu;
6155 matrix[5][6] = sumx1y1sqR;
6156 matrix[5][7] = sumy1cuR;
6157
6158 matrix[6][0] = sumx1R;
6159 matrix[6][1] = sumx1sqR;
6160 matrix[6][2] = sumx1y1R;
6161 matrix[6][3] = sumx1cuR;
6162 matrix[6][4] = sumx1sqy1R;
6163 matrix[6][5] = sumx1y1sqR;
6164 matrix[6][6] = sumx1sqRsq;
6165 matrix[6][7] = sumx1y1Rsq;
6166
6167 matrix[7][0] = sumy1R;
6168 matrix[7][1] = sumx1y1R;
6169 matrix[7][2] = sumy1sqR;
6170 matrix[7][3] = sumx1sqy1R;
6171 matrix[7][4] = sumx1y1sqR;
6172 matrix[7][5] = sumy1cuR;
6173 matrix[7][6] = sumx1y1Rsq;
6174 matrix[7][7] = sumy1sqRsq;
6175
6176 vector[0] = sumy2;
6177 vector[1] = sumy2x1;
6178 vector[2] = sumy2y1;
6179 vector[3] = sumy2x1sq;
6180 vector[4] = sumy2x1y1;
6181 vector[5] = sumy2y1sq;
6182 vector[6] = sumy2x1R;
6183 vector[7] = sumy2y1R;
6184
6185 #ifdef DEBUG
6186 printf("before calling solution routines for IJKLMNOP, here's matrix\n");
6187 print_matrix(matrix, 8);
6188 #endif
6189
6190 /*
6191 * and now call the Gaussian-elimination routines to solve the matrix.
6192 * The solution for TRANS coefficients I, J, K, L, M, N, O, P will be placed
6193 * into the elements on "vector" after "gauss_matrix" finishes.
6194 */
6195 if (gauss_matrix(matrix, 8, vector) != SH_SUCCESS) {
6196 shError("calc_trans_cubic: can't solve for coeffs I,J,K,L,M,N,O,P ");
6197 return(SH_GENERIC_ERROR);
6198 }
6199
6200 #ifdef DEBUG
6201 printf("after calling solution routines, here's matrix\n");
6202 print_matrix(matrix, 8);
6203 #endif
6204
6205 solved_i = vector[0];
6206 solved_j = vector[1];
6207 solved_k = vector[2];
6208 solved_l = vector[3];
6209 solved_m = vector[4];
6210 solved_n = vector[5];
6211 solved_o = vector[6];
6212 solved_p = vector[7];
6213
6214
6215 /*
6216 * assign the coefficients we've just calculated to the output
6217 * TRANS structure.
6218 */
6219 trans->a = solved_a;
6220 trans->b = solved_b;
6221 trans->c = solved_c;
6222 trans->d = solved_d;
6223 trans->e = solved_e;
6224 trans->f = solved_f;
6225 trans->g = solved_g;
6226 trans->h = solved_h;
6227 trans->i = solved_i;
6228 trans->j = solved_j;
6229 trans->k = solved_k;
6230 trans->l = solved_l;
6231 trans->m = solved_m;
6232 trans->n = solved_n;
6233 trans->o = solved_o;
6234 trans->p = solved_p;
6235
6236 /*
6237 * free up memory we allocated for this function
6238 */
6239 free_matrix(matrix, 8);
6240
6241 return(SH_SUCCESS);
6242 }
6243
6244
6245
6246
6247 /***************************************************************************
6248 * PROCEDURE: gauss_matrix
6249 *
6250 * DESCRIPTION:
6251 * Given a square 2-D 'num'-by-'num' matrix, called "matrix", and given
6252 * a 1-D vector "vector" of 'num' elements, find the 1-D vector
6253 * called "solution_vector" which satisfies the equation
6254 *
6255 * matrix * solution_vector = vector
6256 *
6257 * where the * above represents matrix multiplication.
6258 *
6259 * What we do is to use Gaussian elimination (with partial pivoting)
6260 * and back-substitution to find the solution_vector.
6261 * We do not pivot in place, but physically move values -- it
6262 * doesn't take much time in this application. After we have found the
6263 * "solution_vector", we replace the contents of "vector" with the
6264 * "solution_vector".
6265 *
6266 * This is a common algorithm. See any book on linear algebra or
6267 * numerical solutions; for example, "Numerical Methods for Engineers,"
6268 * by Steven C. Chapra and Raymond P. Canale, McGraw-Hill, 1998,
6269 * Chapter 9.
6270 *
6271 * If an error occurs (if the matrix is singular), this prints an error
6272 * message and returns with error code.
6273 *
6274 * RETURN:
6275 * SH_SUCCESS if all goes well
6276 * SH_GENERIC_ERROR if not -- if matrix is singular
6277 *
6278 * </AUTO>
6279 */
6280
6281
6282 static int
gauss_matrix(double ** matrix,int num,double * vector)6283 gauss_matrix
6284 (
6285 double **matrix, /* I/O: the square 2-D matrix we'll invert */
6286 /* will hold inverse matrix on output */
6287 int num, /* I: number of rows and cols in matrix */
6288 double *vector /* I/O: vector which holds "b" values in input */
6289 /* and the solution vector "x" on output */
6290 )
6291 {
6292 int i, j, k;
6293 double *biggest_val;
6294 double *solution_vector;
6295 double factor;
6296 double sum;
6297
6298 #ifdef DEBUG
6299 print_matrix(matrix, num);
6300 #endif
6301
6302 biggest_val = (double *) shMalloc(num*sizeof(double));
6303 solution_vector = (double *) shMalloc(num*sizeof(double));
6304
6305 /*
6306 * step 1: we find the largest value in each row of matrix,
6307 * and store those values in 'biggest_val' array.
6308 * We use this information to pivot the matrix.
6309 */
6310 for (i = 0; i < num; i++) {
6311 biggest_val[i] = fabs(matrix[i][0]);
6312 for (j = 1; j < num; j++) {
6313 if (fabs(matrix[i][j]) > biggest_val[i]) {
6314 biggest_val[i] = fabs(matrix[i][j]);
6315 }
6316 }
6317 if (biggest_val[i] == 0.0) {
6318 shError("gauss_matrix: biggest val in row %d is zero", i);
6319 shFree(biggest_val);
6320 shFree(solution_vector);
6321 return(SH_GENERIC_ERROR);
6322 }
6323 }
6324
6325 /*
6326 * step 2: we use Gaussian elimination to convert the "matrix"
6327 * into a triangular matrix, in which the values of all
6328 * elements below the diagonal is zero.
6329 */
6330 for (i = 0; i < num - 1; i++) {
6331
6332 /* pivot this row (if necessary) */
6333 if (gauss_pivot(matrix, num, vector, biggest_val, i) == SH_GENERIC_ERROR) {
6334 shError("gauss_matrix: singular matrix");
6335 shFree(biggest_val);
6336 shFree(solution_vector);
6337 return(SH_GENERIC_ERROR);
6338 }
6339
6340 if (fabs(matrix[i][i]/biggest_val[i]) < MATRIX_TOL) {
6341 shError("gauss_matrix: Y: row %d has tiny value %f / %f",
6342 i, matrix[i][i], biggest_val[i]);
6343 shFree(biggest_val);
6344 shFree(solution_vector);
6345 return(SH_GENERIC_ERROR);
6346 }
6347
6348 /* we eliminate this variable in all rows below the current one */
6349 for (j = i + 1; j < num; j++) {
6350 factor = matrix[j][i]/matrix[i][i];
6351 for (k = i + 1; k < num; k++) {
6352 matrix[j][k] -= factor*matrix[i][k];
6353 }
6354 /* and in the vector, too */
6355 vector[j] -= factor*vector[i];
6356 }
6357
6358 }
6359
6360 /*
6361 * make sure that the last row's single remaining element
6362 * isn't too tiny
6363 */
6364 if (fabs(matrix[num-1][num-1]/biggest_val[num-1]) < MATRIX_TOL) {
6365 shError("gauss_matrix: Z: row %d has tiny value %f / %f",
6366 num, matrix[num-1][num-1], biggest_val[num-1]);
6367 shFree(biggest_val);
6368 shFree(solution_vector);
6369 return(SH_GENERIC_ERROR);
6370 }
6371
6372 /*
6373 * step 3: we can now calculate the solution_vector values
6374 * via back-substitution; we start at the last value in the
6375 * vector (at the "bottom" of the vector) and work
6376 * upwards towards the top.
6377 */
6378 solution_vector[num-1] = vector[num-1] / matrix[num-1][num-1];
6379 for (i = num - 2; i >= 0; i--) {
6380 sum = 0.0;
6381 for (j = i + 1; j < num; j++) {
6382 sum += matrix[i][j]*solution_vector[j];
6383 }
6384 solution_vector[i] = (vector[i] - sum) / matrix[i][i];
6385 }
6386
6387
6388 /*
6389 * step 4: okay, we've found the values in the solution vector!
6390 * We now replace the input values in 'vector' with these
6391 * solution_vector values, and we're done.
6392 */
6393 for (i = 0; i < num; i++) {
6394 vector[i] = solution_vector[i];
6395 }
6396
6397
6398 /* clean up */
6399 shFree(solution_vector);
6400 shFree(biggest_val);
6401
6402 return(SH_SUCCESS);
6403 }
6404
6405
6406 /***************************************************************************
6407 * PROCEDURE: gauss_pivot
6408 *
6409 * DESCRIPTION:
6410 * This routine is called by "gauss_matrix". Given a square "matrix"
6411 * of "num"-by-"num" elements, and given a "vector" of "num" elements,
6412 * and given a particular "row" value, this routine finds the largest
6413 * value in the matrix at/below the given "row" position. If that
6414 * largest value isn't in the given "row", this routine switches
6415 * rows in the matrix (and in the vector) so that the largest value
6416 * will now be in "row".
6417 *
6418 * RETURN:
6419 * SH_SUCCESS if all goes well
6420 * SH_GENERIC_ERROR if not -- if matrix is singular
6421 *
6422 * </AUTO>
6423 */
6424
6425 #define SWAP(a,b) { double temp = (a); (a) = (b); (b) = temp; }
6426
6427 static int
gauss_pivot(double ** matrix,int num,double * vector,double * biggest_val,int row)6428 gauss_pivot
6429 (
6430 double **matrix, /* I/O: a square 2-D matrix we are inverting */
6431 int num, /* I: number of rows and cols in matrix */
6432 double *vector, /* I/O: vector which holds "b" values in input */
6433 double *biggest_val, /* I: largest value in each row of matrix */
6434 int row /* I: want to pivot around this row */
6435 )
6436 {
6437 int i, j;
6438 int col, pivot_row;
6439 double big, other_big;
6440
6441 /* sanity checks */
6442 shAssert(matrix != NULL);
6443 shAssert(vector != NULL);
6444 shAssert(row < num);
6445
6446
6447 pivot_row = row;
6448 big = fabs(matrix[row][row]/biggest_val[row]);
6449 #ifdef DEBUG
6450 print_matrix(matrix, num);
6451 printf(" biggest_val: ");
6452 for (i = 0; i < num; i++) {
6453 printf("%9.4e ", biggest_val[i]);
6454 }
6455 printf("\n");
6456 printf(" gauss_pivot: row %3d %9.4e %9.4e %12.5e ",
6457 row, matrix[row][row], biggest_val[row], big);
6458 #endif
6459
6460 for (i = row + 1; i < num; i++) {
6461 other_big = fabs(matrix[i][row]/biggest_val[i]);
6462 if (other_big > big) {
6463 big = other_big;
6464 pivot_row = i;
6465 }
6466 }
6467 #ifdef DEBUG
6468 printf(" pivot_row %3d %9.4e %9.4e %12.5e ",
6469 pivot_row, matrix[pivot_row][pivot_row],
6470 biggest_val[pivot_row], big);
6471 #endif
6472
6473 /*
6474 * if another row is better for pivoting, switch it with 'row'
6475 * and switch the corresponding elements in 'vector'
6476 * and switch the corresponding elements in 'biggest_val'
6477 */
6478 if (pivot_row != row) {
6479 #ifdef DEBUG
6480 printf(" will swap \n");
6481 #endif
6482 for (col = row; col < num; col++) {
6483 SWAP(matrix[pivot_row][col], matrix[row][col]);
6484 }
6485 SWAP(vector[pivot_row], vector[row]);
6486 SWAP(biggest_val[pivot_row], biggest_val[row]);
6487 }
6488 else {
6489 #ifdef DEBUG
6490 printf(" no swap \n");
6491 #endif
6492 }
6493
6494 return(SH_SUCCESS);
6495 }
6496
6497
6498
6499 /***************************************************************************
6500 * PROCEDURE: atFindMedtf
6501 *
6502 * DESCRIPTION:
6503 * Assume that the two input lists of stars were taken by the same
6504 * instrument, so that they have the same scale and rotation.
6505 * If this is true, then a simple translation ought to register
6506 * the two lists. This routine tries to characterize that
6507 * translation: it finds both the mean shift in (x, y),
6508 * and the median shift in (x, y). It also calculates the
6509 * (clipped) standard deviation from the mean shift. The results of all
6510 * these calculations are placed into the given MEDTF structure.
6511 *
6512 * RETURN:
6513 * SH_SUCCESS if all goes well
6514 * SH_GENERIC_ERROR if not
6515 *
6516 * </AUTO>
6517 */
6518
6519 int
atFindMedtf(int num_matched_A,struct s_star * listA,int num_matched_B,struct s_star * listB,double medsigclip,MEDTF * medtf)6520 atFindMedtf
6521 (
6522 int num_matched_A, /* I: number of matched stars in list A */
6523 struct s_star *listA, /* I: list of matched stars from set A */
6524 int num_matched_B, /* I: number of matched stars in list B */
6525 struct s_star *listB, /* I: list of matched stars from set B */
6526 double medsigclip, /* I: sigma-clipping factor */
6527 MEDTF *medtf /* O: we place results into this structure */
6528 )
6529 {
6530 int i, nstar, num_within_clip;
6531 double *dx, *dy, *dxclip, *dyclip;
6532 double xdist, ydist;
6533 double Dx_sum, Dx_sum2, Dx_rms, Dx_ave, Dx_med;
6534 double Dy_sum, Dy_sum2, Dy_rms, Dy_ave, Dy_med, clip;
6535 struct s_star *A_star, *B_star;
6536
6537
6538
6539 if (num_matched_A < 3) {
6540 shError("atFindMedtf: fewer than 3 matched pairs; cannot find MEDTF");
6541 return(SH_GENERIC_ERROR);
6542 }
6543 nstar = num_matched_A;
6544
6545 /* sanity checks */
6546 shAssert(num_matched_A == num_matched_B);
6547 shAssert(listA != NULL);
6548 shAssert(listB != NULL);
6549 shAssert(medsigclip >= 0);
6550 shAssert(medtf != NULL);
6551
6552 /*
6553 * allocate space for arrays we use to sort distances
6554 * between matched items in the lists
6555 */
6556 dx = (double *) shMalloc(nstar*sizeof(double));
6557 dy = (double *) shMalloc(nstar*sizeof(double));
6558
6559 /*
6560 * Step 1: calculate distances between matched stars,
6561 * and fill the "dx" and "dy" arrays for later calculation
6562 * of the median
6563 */
6564 Dx_sum = 0.0;
6565 Dy_sum = 0.0;
6566 Dx_sum2 = 0.0;
6567 Dy_sum2 = 0.0;
6568 A_star = listA;
6569 B_star = listB;
6570 for (i = 0; i < nstar; i++) {
6571 shAssert(A_star != NULL);
6572 shAssert(B_star != NULL);
6573
6574 xdist = (double) B_star->x - (double) A_star->x;
6575 ydist = (double) B_star->y - (double) A_star->y;
6576 dx[i] = xdist;
6577 dy[i] = ydist;
6578 #ifdef DEBUG
6579 printf (" medtf: %4d %4d xa %7.4f xb %7.4f ya %7.4f yb %7.4f\n",
6580 A_star->id, B_star->id,
6581 (double) A_star->x, (double) B_star->x,
6582 (double) A_star->y, (double) B_star->y);
6583 printf (" medtf: %4d dx %7.4f dy %7.4f \n", i, xdist, ydist);
6584 #endif
6585 Dx_sum += xdist;
6586 Dy_sum += ydist;
6587 Dx_sum2 += xdist*xdist;
6588 Dy_sum2 += ydist*ydist;
6589
6590 A_star = A_star->next;
6591 B_star = B_star->next;
6592 }
6593
6594 /*
6595 * Step 2: calculate the mean distances and (unclipped) stdev
6596 */
6597 Dx_ave = Dx_sum/nstar;
6598 Dy_ave = Dy_sum/nstar;
6599 Dx_rms = sqrt(Dx_sum2/nstar - Dx_ave*Dx_ave);
6600 Dy_rms = sqrt(Dy_sum2/nstar - Dy_ave*Dy_ave);
6601
6602 /*
6603 * Step 3: calculate the median distances
6604 */
6605 qsort((char *) dx, nstar, sizeof(double), (PFI) compare_double);
6606 Dx_med = find_percentile(dx, nstar, (double) 0.50);
6607 qsort((char *) dy, nstar, sizeof(double), (PFI) compare_double);
6608 Dy_med = find_percentile(dy, nstar, (double) 0.50);
6609
6610 /*
6611 * Step 4 (if desired): recalculate statistics, using only pairs of
6612 * stars separated by 'medsigclip' standard
6613 * deviations from the mean.
6614 */
6615 if (medsigclip > 0) {
6616 if ((Dx_rms <= 0.0) || (Dy_rms <= 0.0)) {
6617 shError("atFindMedtf: RMS <= 0.0, so can't calculate clipped values");
6618 }
6619 else {
6620
6621 /*
6622 * we need another pair of arrays, this time containing only
6623 * those distances within the clipping criterion
6624 */
6625 dxclip = (double *) shMalloc(nstar*sizeof(double));
6626 dyclip = (double *) shMalloc(nstar*sizeof(double));
6627
6628 /* calculate the maximum acceptable distance */
6629 clip = medsigclip*sqrt(Dx_rms*Dy_rms);
6630 if ((fabs(Dx_med - Dx_ave) > 0.5*clip) ||
6631 (fabs(Dy_med - Dy_ave) > 0.5*clip)) {
6632 shError("atFindMedtf: dangerous skewness in shifts");
6633 }
6634
6635 /*
6636 * recalculate statistics, discarding values outside the
6637 * clipping criterion
6638 */
6639 Dx_sum = 0.0;
6640 Dy_sum = 0.0;
6641 Dx_sum2 = 0.0;
6642 Dy_sum2 = 0.0;
6643 num_within_clip = 0;
6644 for (i = 0; i < nstar; i++) {
6645 if (fabs(dx[i] - Dx_med) > clip) {
6646 continue;
6647 }
6648 if (fabs(dy[i] - Dy_med) > clip) {
6649 continue;
6650 }
6651
6652 xdist = dx[i];
6653 ydist = dy[i];
6654 dxclip[num_within_clip] = xdist;
6655 dyclip[num_within_clip] = ydist;
6656 Dx_sum += xdist;
6657 Dy_sum += ydist;
6658 Dx_sum2 += xdist*xdist;
6659 Dy_sum2 += ydist*ydist;
6660 num_within_clip++;
6661 }
6662
6663 Dx_ave = Dx_sum/num_within_clip;
6664 Dy_ave = Dy_sum/num_within_clip;
6665 Dx_rms = sqrt(Dx_sum2/num_within_clip - Dx_ave*Dx_ave);
6666 Dy_rms = sqrt(Dy_sum2/num_within_clip - Dy_ave*Dy_ave);
6667
6668 qsort((char *) dxclip, num_within_clip, sizeof(double),
6669 (PFI) compare_double);
6670 Dx_med = find_percentile(dxclip, num_within_clip, (double) 0.50);
6671 qsort((char *) dyclip, num_within_clip, sizeof(double),
6672 (PFI) compare_double);
6673 Dy_med = find_percentile(dyclip, num_within_clip, (double) 0.50);
6674
6675 shFree(dxclip);
6676 shFree(dyclip);
6677
6678 nstar = num_within_clip;
6679 }
6680 }
6681
6682 /* finally, set the values in the MEDTF structure */
6683 medtf->mdx = Dx_med;
6684 medtf->mdy = Dy_med;
6685 medtf->adx = Dx_ave;
6686 medtf->ady = Dy_ave;
6687 medtf->sdx = Dx_rms;
6688 medtf->sdy = Dy_rms;
6689 medtf->nm = nstar;
6690
6691 /* free the memory we've allocated */
6692 shFree(dx);
6693 shFree(dy);
6694
6695 return(SH_SUCCESS);
6696 }
6697
6698
6699 /***************************************************************************
6700 * PROCEDURE: atCalcRMS
6701 * Added 17 Jan 2002 by JPB
6702 *
6703 * DESCRIPTION:
6704 * This routine takes two matched lists and calculates the
6705 * rms of the differences. Just for simple diagnostic purposes.
6706 *
6707 * If the two lists are empty, set the output args to 0.00 and
6708 * return with code SH_SUCCESS.
6709 *
6710 * RETURN:
6711 * SH_SUCCESS if all goes well (even if lists were empty)
6712 * SH_GENERIC_ERROR if not
6713 *
6714 */
atCalcRMS(int num_A,struct s_star * mlistA,int num_B,struct s_star * mlistB,double * Dx_rms,double * Dy_rms)6715 int atCalcRMS
6716 (
6717 int num_A, /* I: number of matched stars in list A */
6718 struct s_star *mlistA, /* I: list of matched stars from set A */
6719 int num_B, /* I: number of matched stars in list B */
6720 struct s_star *mlistB, /* I: list of matched stars from set B */
6721 double *Dx_rms,
6722 double *Dy_rms
6723 )
6724 {
6725 double Dxterm, Dyterm, Dx_sum2, Dy_sum2, xms, yms;
6726 int ii, Nstar, Ntoss;
6727 struct s_star *A_star, *B_star;
6728
6729 shAssert(num_A == num_B);
6730 if (num_A == 0) {
6731 *Dx_rms = 0.0;
6732 *Dy_rms = 0.0;
6733 return(SH_SUCCESS);
6734 }
6735
6736 shAssert(mlistA != NULL);
6737 shAssert(mlistB != NULL);
6738 Nstar = num_A;
6739 Dx_sum2 = Dy_sum2 = 0.0;
6740
6741 for (ii = 0, A_star = mlistA, B_star = mlistB; ii < Nstar;
6742 ii++, A_star = A_star->next, B_star = B_star->next) {
6743 shAssert(A_star != NULL);
6744 shAssert(B_star != NULL);
6745 Dxterm = (double)(B_star->x - A_star->x);
6746 Dyterm = (double)(B_star->y - A_star->y);
6747 Dx_sum2 += Dxterm*Dxterm;
6748 Dy_sum2 += Dyterm*Dyterm;
6749 }
6750 xms = (Dx_sum2/Nstar); /* these are mean-squared! */
6751 yms = (Dy_sum2/Nstar);
6752
6753 /* Ok, we just do a quick conservative 3-sigma clip here */
6754
6755 Dx_sum2 = Dy_sum2 = 0.0;
6756 for(Ntoss = ii = 0, A_star = mlistA, B_star = mlistB; ii < Nstar;
6757 ii++, A_star = A_star->next, B_star = B_star->next) {
6758 Dxterm = (double)(B_star->x - A_star->x);
6759 Dyterm = (double)(B_star->y - A_star->y);
6760
6761 /* Note: squaring these terms, unlike loop above! */
6762 Dxterm *= Dxterm;
6763 Dyterm *= Dyterm;
6764
6765 if (Dxterm < 9*xms && Dyterm < 9*yms) {
6766 Dx_sum2 += Dxterm;
6767 Dy_sum2 += Dyterm;
6768 }
6769 else {
6770 Ntoss++;
6771 }
6772 }
6773 if (Dx_sum2 <= 0.0) {
6774 *Dx_rms = 0.0;
6775 }
6776 else {
6777 *Dx_rms = sqrt(Dx_sum2/(Nstar-Ntoss));
6778 }
6779 if (Dy_sum2 <= 0.0) {
6780 *Dy_rms = 0.0;
6781 }
6782 else {
6783 *Dy_rms = sqrt(Dy_sum2/(Nstar-Ntoss));
6784 }
6785
6786 return(SH_SUCCESS);
6787 }
6788
6789
6790
6791
6792 /***************************************************************************
6793 * PROCEDURE: is_desired_rotation
6794 *
6795 * DESCRIPTION:
6796 * We are given two triangles and a desired relative rotation angle
6797 * (and tolerance).
6798 * Our job is to determine if the two triangles are rotated relative
6799 * to each other by the given angle, within the tolerance.
6800 *
6801 * This is a little tricky, because we need to allow for angles
6802 * wrapping around the range of (-pi, +pi) in radians, or (-180, +180)
6803 * in degrees. We therefore make two checks, first without any
6804 * wrapping, and then with the appropriate wrapping; if either check
6805 * succeeds, we return with SH_SUCCESS.
6806 *
6807 * We place the actual measured rotation angle (in degrees) between
6808 * these two triangles into the 'actual_angle' argument.
6809 *
6810 * RETURN:
6811 * 1 if the two triangles are rotated by the desired angle
6812 * 0 if not
6813 *
6814 * </AUTO>
6815 */
6816
6817 static int
is_desired_rotation(struct s_triangle * tri_A,struct s_triangle * tri_B,double want_angle_deg,double tolerance_deg,double * actual_angle_deg)6818 is_desired_rotation
6819 (
6820 struct s_triangle *tri_A, /* I: first triangle */
6821 struct s_triangle *tri_B, /* I: second triangle */
6822 double want_angle_deg, /* I: desired rotation, in degrees */
6823 double tolerance_deg, /* I: accept any rotation within this amount */
6824 /* I: of desired value, in degrees */
6825 double *actual_angle_deg /* O: the measured rotation angle between */
6826 /* the triangles, in degrees */
6827 )
6828 {
6829 int is_good_angle;
6830 double min_angle_deg, max_angle_deg, wrapped_delta_deg;
6831 double delta_angle, delta_angle_deg;
6832
6833 is_good_angle = 0;
6834 min_angle_deg = want_angle_deg - tolerance_deg;
6835 max_angle_deg = want_angle_deg + tolerance_deg;
6836
6837 delta_angle = tri_A->side_a_angle - tri_B->side_a_angle;
6838 delta_angle_deg = delta_angle*(180.0/3.15159);
6839 #ifdef DEBUG2
6840 printf("delta_angle_deg %7.2f ", delta_angle_deg);
6841 #endif
6842
6843 if ((delta_angle_deg >= min_angle_deg) && (delta_angle_deg <= max_angle_deg)) {
6844 #ifdef DEBUG2
6845 printf(" is good ");
6846 #endif
6847 is_good_angle = 1;
6848 }
6849 #ifdef DEBUG2
6850 else {
6851 printf(" is bad ");
6852 }
6853 #endif
6854
6855 /*
6856 * now we check to see if the wrapped version of the angle falls
6857 * into the desired range
6858 */
6859 if (delta_angle_deg > 0) {
6860 wrapped_delta_deg = delta_angle_deg - 360.0;
6861 }
6862 else {
6863 wrapped_delta_deg = delta_angle_deg + 360.0;
6864 }
6865 #ifdef DEBUG2
6866 printf("wrapped_delta_deg %7.2f ", wrapped_delta_deg);
6867 #endif
6868 if ((wrapped_delta_deg >= min_angle_deg) && (wrapped_delta_deg <= max_angle_deg)) {
6869 #ifdef DEBUG2
6870 printf(" is good\n");
6871 #endif
6872 is_good_angle = 1;
6873 }
6874 #ifdef DEBUG2
6875 else {
6876 printf(" is bad \n");
6877 }
6878 #endif
6879
6880
6881 *actual_angle_deg = delta_angle_deg;
6882
6883 if (is_good_angle == 1) {
6884 return(1);
6885 }
6886 else {
6887 return(0);
6888 }
6889 }
6890
6891