1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5
6 //# Modified by Robert Lancaster for the StellarSolver Internal Library
7 #ifdef _MSC_VER
8 #define _USE_MATH_DEFINES
9 #endif
10
11 #include <stdio.h>
12 #include <string.h>
13 #include <math.h>
14 #include <assert.h>
15 #include <sys/types.h>
16 #include <unistd.h>
17 #include <stdarg.h>
18
19 #include "os-features.h"
20 #include "ioutils.h"
21 #include "mathutil.h"
22 #include "matchobj.h"
23 #include "solver.h"
24 #include "verify.h"
25 #include "tic.h"
26 #include "solvedclient.h"
27 #include "solvedfile.h"
28 #include "fit-wcs.h"
29 #include "sip-utils.h"
30 #include "keywords.h"
31 #include "log.h"
32 #include "pquad.h"
33 #include "kdtree.h"
34 #include "quad-utils.h"
35 #include "errors.h"
36 #include "tweak2.h"
37
38 #if TESTING_TRYALLCODES
39 #define DEBUGSOLVER 1
40 #define TRY_ALL_CODES test_try_all_codes
41 void test_try_all_codes(pquad* pq,
42 int* fieldstars, int dimquad,
43 solver_t* solver, double tol2);
44
45 #else
46 #define TRY_ALL_CODES try_all_codes
47 #endif
48
49 #if TESTING_TRYPERMUTATIONS
50 #define DEBUGSOLVER 1
51 #define TEST_TRY_PERMUTATIONS test_try_permutations
52 void test_try_permutations(int* stars, double* code, int dimquad, solver_t* s);
53
54 #else
55 #define TEST_TRY_PERMUTATIONS(u,v,x,y) // no-op.
56 #endif
57
58
59
60
61
62
solver_set_keep_logodds(solver_t * solver,double logodds)63 void solver_set_keep_logodds(solver_t* solver, double logodds) {
64 solver->logratio_tokeep = logodds;
65 }
66
solver_set_parity(solver_t * solver,int parity)67 int solver_set_parity(solver_t* solver, int parity) {
68 if (!((parity == PARITY_NORMAL) || (parity == PARITY_FLIP) || (parity == PARITY_BOTH))) {
69 ERROR("Invalid parity value: %i", parity);
70 return -1;
71 }
72 solver->parity = parity;
73 return 0;
74 }
75
solver_did_solve(const solver_t * solver)76 anbool solver_did_solve(const solver_t* solver) {
77 return solver->best_match_solves;
78 }
79
solver_get_quad_size_range_arcsec(const solver_t * solver,double * qmin,double * qmax)80 void solver_get_quad_size_range_arcsec(const solver_t* solver, double* qmin, double* qmax) {
81 if (qmin) {
82 *qmin = solver->quadsize_min * solver_get_pixscale_low(solver);
83 }
84 if (qmax) {
85 double q = solver->quadsize_max;
86 if (q == 0)
87 q = solver->field_diag;
88 *qmax = q * solver_get_pixscale_high(solver);
89 }
90 }
91
solver_get_field_jitter(const solver_t * solver)92 double solver_get_field_jitter(const solver_t* solver) {
93 return solver->verify_pix;
94 }
95
solver_get_field_center(const solver_t * solver,double * px,double * py)96 void solver_get_field_center(const solver_t* solver, double* px, double* py) {
97 if (px)
98 *px = (solver->field_maxx + solver->field_minx)/2.0;
99 if (py)
100 *py = (solver->field_maxy + solver->field_miny)/2.0;
101 }
102
solver_get_max_radius_arcsec(const solver_t * solver)103 double solver_get_max_radius_arcsec(const solver_t* solver) {
104 return solver->funits_upper * solver->field_diag / 2.0;
105 }
106
solver_get_best_match(solver_t * solver)107 MatchObj* solver_get_best_match(solver_t* solver) {
108 return &(solver->best_match);
109 }
110
solver_get_best_match_index_name(const solver_t * solver)111 const char* solver_get_best_match_index_name(const solver_t* solver) {
112 return solver->best_index->indexname;
113 }
114
solver_get_pixscale_low(const solver_t * solver)115 double solver_get_pixscale_low(const solver_t* solver) {
116 return solver->funits_lower;
117 }
solver_get_pixscale_high(const solver_t * solver)118 double solver_get_pixscale_high(const solver_t* solver) {
119 return solver->funits_upper;
120 }
121
solver_set_quad_size_range(solver_t * solver,double qmin,double qmax)122 void solver_set_quad_size_range(solver_t* solver, double qmin, double qmax) {
123 solver->quadsize_min = qmin;
124 solver->quadsize_max = qmax;
125 }
126
solver_set_quad_size_fraction(solver_t * solver,double qmin,double qmax)127 void solver_set_quad_size_fraction(solver_t* solver, double qmin, double qmax) {
128 solver_set_quad_size_range(solver, qmin * MIN(solver_field_width(solver), solver_field_height(solver)),
129 qmax * solver->field_diag);
130 }
131
solver_tweak2(solver_t * sp,MatchObj * mo,int order,sip_t * verifysip)132 void solver_tweak2(solver_t* sp, MatchObj* mo, int order, sip_t* verifysip) {
133 double* xy = NULL;
134 int Nxy;
135 double indexjitter;
136 // quad center
137 double qc[2];
138 // quad radius-squared
139 double Q2;
140 // initial WCS
141 sip_t startsip;
142 int* theta;
143 double* odds;
144 double* refradec;
145 int i;
146 double newodds;
147 int nm, nc, nd;
148 int besti;
149 int startorder;
150
151 indexjitter = mo->index_jitter; // ref cat positional error, in arcsec.
152 xy = starxy_to_xy_array(sp->fieldxy, NULL);
153 Nxy = starxy_n(sp->fieldxy);
154 qc[0] = (mo->quadpix[0] + mo->quadpix[2]) / 2.0;
155 qc[1] = (mo->quadpix[1] + mo->quadpix[3]) / 2.0;
156 Q2 = 0.25 * distsq(mo->quadpix, mo->quadpix + 2, 2);
157 if (Q2 == 0.0) {
158 // can happen if we're verifying an existing WCS
159 // note, this is radius-squared, so 1e6 is not crazy.
160 Q2 = 1e6;
161 // set qc to the image center here? or crpix?
162 logverb("solver_tweak2(): setting Q2=%g; qc=(%g,%g)\n", Q2, qc[0], qc[1]);
163 }
164
165 // mo->refradec may be NULL at this point, so get it from refxyz instead...
166 refradec = malloc(3 * mo->nindex * sizeof(double));
167 for (i=0; i<mo->nindex; i++)
168 xyzarr2radecdegarr(mo->refxyz + i*3, refradec + i*2);
169
170 // Verifying an existing WCS?
171 if (verifysip) {
172 memcpy(&startsip, verifysip, sizeof(sip_t));
173 startorder = MIN(verifysip->a_order, sp->tweak_aborder);
174 } else {
175 startorder = 1;
176 sip_wrap_tan(&(mo->wcstan), &startsip);
177 }
178
179 startsip.ap_order = startsip.bp_order = sp->tweak_abporder;
180 startsip.a_order = startsip.b_order = sp->tweak_aborder;
181 logverb("solver_tweak2: setting orders %i, %i\n", sp->tweak_aborder, sp->tweak_abporder);
182
183 // for TWEAK_DEBUG_PLOTs
184 theta = mo->theta;
185 besti = mo->nbest-1;//mo->nmatch + mo->nconflict + mo->ndistractor;
186
187 logverb("solver_tweak2: set_crpix %i, crpix (%.1f,%.1f)\n",
188 sp->set_crpix, sp->crpix[0], sp->crpix[1]);
189 mo->sip = tweak2(xy, Nxy,
190 sp->verify_pix, // pixel positional noise sigma
191 solver_field_width(sp),
192 solver_field_height(sp),
193 refradec, mo->nindex,
194 indexjitter, qc, Q2,
195 sp->distractor_ratio,
196 sp->logratio_bail_threshold,
197 order, sp->tweak_abporder,
198 &startsip, NULL, &theta, &odds,
199 sp->set_crpix ? sp->crpix : NULL,
200 &newodds, &besti, mo->testperm, startorder);
201 free(refradec);
202
203 // FIXME -- update refxy? Nobody uses it, right?
204 free(mo->refxy);
205 mo->refxy = NULL;
206 // FIXME -- and testperm?
207 free(mo->testperm);
208 mo->testperm = NULL;
209
210 if (mo->sip) {
211 // Yoink the TAN solution (?)
212 memcpy(&(mo->wcstan), &(mo->sip->wcstan), sizeof(tan_t));
213
214 // Plug in the new "theta" and "odds".
215 free(mo->theta);
216 free(mo->matchodds);
217 mo->theta = theta;
218 mo->matchodds = odds;
219
220 mo->logodds = newodds;
221
222 verify_count_hits(theta, besti, &nm, &nc, &nd);
223 mo->nmatch = nm;
224 mo->nconflict = nc;
225 mo->ndistractor = nd;
226 matchobj_compute_derived(mo);
227 }
228 free(xy);
229 }
230
solver_log_params(const solver_t * sp)231 void solver_log_params(const solver_t* sp) {
232 int i;
233 logverb("Solver:\n");
234 logverb(" Arcsec per pix range: %g, %g\n", sp->funits_lower, sp->funits_upper);
235 logverb(" Image size: %g x %g\n", solver_field_width(sp), solver_field_height(sp));
236 logverb(" Quad size range: %g, %g\n", sp->quadsize_min, sp->quadsize_max);
237 logverb(" Objs: %i, %i\n", sp->startobj, sp->endobj);
238 logverb(" Parity: %i, %s\n", sp->parity, sp->parity == PARITY_NORMAL ? "normal" : (sp->parity == PARITY_FLIP ? "flip" : "both"));
239 if (sp->use_radec) {
240 double ra,dec,rad;
241 xyzarr2radecdeg(sp->centerxyz, &ra, &dec);
242 rad = distsq2deg(sp->r2);
243 logverb(" Use_radec? yes, (%g, %g), radius %g deg\n", ra, dec, rad);
244 } else {
245 logverb(" Use_radec? no\n");
246 }
247 logverb(" Verify_pix: %g\n", sp->verify_pix);
248 logverb(" Code tol: %g\n", sp->codetol);
249 logverb(" Dist from quad bonus: %s\n", sp->distance_from_quad_bonus ? "yes" : "no");
250 logverb(" Distractor ratio: %g\n", sp->distractor_ratio);
251 logverb(" Log tune-up threshold: %g\n", sp->logratio_totune);
252 logverb(" Log bail threshold: %g\n", sp->logratio_bail_threshold);
253 logverb(" Log stoplooking threshold: %g\n", sp->logratio_stoplooking);
254 logverb(" Maxquads %i\n", sp->maxquads);
255 logverb(" Maxmatches %i\n", sp->maxmatches);
256 logverb(" Set CRPIX? %s", sp->set_crpix ? "yes" : "no\n");
257 if (sp->set_crpix) {
258 if (sp->set_crpix_center)
259 logverb(", center\n");
260 else
261 logverb(", %g, %g\n", sp->crpix[0], sp->crpix[1]);
262 }
263 logverb(" Tweak? %s\n", sp->do_tweak ? "yes" : "no");
264 if (sp->do_tweak) {
265 logverb(" Forward order %i\n", sp->tweak_aborder);
266 logverb(" Reverse order %i\n", sp->tweak_abporder);
267 }
268 logverb(" Indexes: %zu\n", pl_size(sp->indexes));
269 for (i=0; i<pl_size(sp->indexes); i++) {
270 index_t* ind = pl_get(sp->indexes, i);
271 logverb(" %s\n", ind->indexname);
272 }
273 logverb(" Field: %i stars\n", starxy_n(sp->fieldxy));
274 for (i=0; i<starxy_n(sp->fieldxy); i++) {
275 debug(" xy (%.1f, %.1f), flux %.1f\n",
276 starxy_getx(sp->fieldxy, i), starxy_gety(sp->fieldxy, i),
277 sp->fieldxy->flux ? starxy_get_flux(sp->fieldxy, i) : 0.0);
278 }
279 }
280
281
solver_print_to(const solver_t * sp,FILE * stream)282 void solver_print_to(const solver_t* sp, FILE* stream) {
283 //int oldlevel = log_get_level();
284 FILE* oldfid = log_get_fid();
285 //log_set_level(LOG_ALL);
286 log_to(stream);
287 solver_log_params(sp);
288 //log_set_level(oldlevel);
289 log_to(oldfid);
290 }
291
292 /*
293 static MatchObj* matchobj_copy_deep(const MatchObj* mo, MatchObj* dest) {
294 if (!dest)
295 dest = malloc(sizeof(MatchObj));
296 memcpy(dest, mo, sizeof(MatchObj));
297 // various modules add things to a mo...
298 blind_matchobj_deep_copy(mo, dest);
299 verify_matchobj_deep_copy(mo, dest);
300 return dest;
301 }
302
303 static void matchobj_free_data(MatchObj* mo) {
304 verify_free_matchobj(mo);
305 blind_free_matchobj(mo);
306 }
307 */
308
309 static const int A = 0, B = 1, C = 2, D = 3;
310
311 // Number of stars in the "backbone" of the quad: stars A and B.
312 static const int NBACK = 2;
313
314 static void find_field_boundaries(solver_t* solver);
315
getx(const double * d,int ind)316 static inline double getx(const double* d, int ind) {
317 return d[ind*2];
318 }
gety(const double * d,int ind)319 static inline double gety(const double* d, int ind) {
320 return d[ind*2 + 1];
321 }
setx(double * d,int ind,double val)322 static inline void setx(double* d, int ind, double val) {
323 d[ind*2] = val;
324 }
sety(double * d,int ind,double val)325 static inline void sety(double* d, int ind, double val) {
326 d[ind*2 + 1] = val;
327 }
328
field_getxy(solver_t * sp,int index,double * x,double * y)329 static void field_getxy(solver_t* sp, int index, double* x, double* y) {
330 *x = starxy_getx(sp->fieldxy, index);
331 *y = starxy_gety(sp->fieldxy, index);
332 }
333
field_getx(solver_t * sp,int index)334 static double field_getx(solver_t* sp, int index) {
335 return starxy_getx(sp->fieldxy, index);
336 }
field_gety(solver_t * sp,int index)337 static double field_gety(solver_t* sp, int index) {
338 return starxy_gety(sp->fieldxy, index);
339 }
340
update_timeused(solver_t * sp)341 static void update_timeused(solver_t* sp) {
342 double usertime, systime;
343 get_resource_stats(&usertime, &systime, NULL);
344 sp->timeused = (usertime + systime) - sp->starttime;
345 if (sp->timeused < 0.0)
346 sp->timeused = 0.0;
347 }
348
set_matchobj_template(solver_t * solver,MatchObj * mo)349 static void set_matchobj_template(solver_t* solver, MatchObj* mo) {
350 if (solver->mo_template)
351 memcpy(mo, solver->mo_template, sizeof(MatchObj));
352 else
353 memset(mo, 0, sizeof(MatchObj));
354 }
355
get_field_center(solver_t * s,double * cx,double * cy)356 static void get_field_center(solver_t* s, double* cx, double* cy) {
357 *cx = 0.5 * (s->field_minx + s->field_maxx);
358 *cy = 0.5 * (s->field_miny + s->field_maxy);
359 }
360
get_field_ll_corner(solver_t * s,double * lx,double * ly)361 static void get_field_ll_corner(solver_t* s, double* lx, double* ly) {
362 *lx = s->field_minx;
363 *ly = s->field_miny;
364 }
365
solver_reset_counters(solver_t * s)366 void solver_reset_counters(solver_t* s) {
367 s->quit_now = FALSE;
368 s->have_best_match = FALSE;
369 s->best_match_solves = FALSE;
370 s->numtries = 0;
371 s->nummatches = 0;
372 s->numscaleok = 0;
373 s->last_examined_object = 0;
374 s->num_cxdx_skipped = 0;
375 s->num_radec_skipped = 0;
376 s->num_abscale_skipped = 0;
377 s->num_verified = 0;
378 }
379
solver_field_width(const solver_t * s)380 double solver_field_width(const solver_t* s) {
381 return s->field_maxx - s->field_minx;
382 }
solver_field_height(const solver_t * s)383 double solver_field_height(const solver_t* s) {
384 return s->field_maxy - s->field_miny;
385 }
386
solver_set_radec(solver_t * s,double ra,double dec,double radius_deg)387 void solver_set_radec(solver_t* s, double ra, double dec, double radius_deg) {
388 s->use_radec = TRUE;
389 radecdeg2xyzarr(ra, dec, s->centerxyz);
390 s->r2 = deg2distsq(radius_deg);
391 }
392
solver_clear_radec(solver_t * s)393 void solver_clear_radec(solver_t* s) {
394 s->use_radec = FALSE;
395 }
396
set_center_and_radius(solver_t * solver,MatchObj * mo,tan_t * tan,sip_t * sip)397 static void set_center_and_radius(solver_t* solver, MatchObj* mo,
398 tan_t* tan, sip_t* sip) {
399 double cx, cy, lx, ly;
400 double xyz[3];
401 get_field_center(solver, &cx, &cy);
402 get_field_ll_corner(solver, &lx, &ly);
403 if (sip) {
404 sip_pixelxy2xyzarr(sip, cx, cy, mo->center);
405 sip_pixelxy2xyzarr(sip, lx, ly, xyz);
406 } else {
407 tan_pixelxy2xyzarr(tan, cx, cy, mo->center);
408 tan_pixelxy2xyzarr(tan, lx, ly, xyz);
409 }
410 mo->radius = sqrt(distsq(mo->center, xyz, 3));
411 mo->radius_deg = dist2deg(mo->radius);
412 }
413
set_index(solver_t * s,index_t * index)414 static void set_index(solver_t* s, index_t* index) {
415 s->index = index;
416 s->rel_index_noise2 = square(index->index_jitter / index->index_scale_lower);
417 }
418
set_diag(solver_t * s)419 static void set_diag(solver_t* s) {
420 s->field_diag = hypot(solver_field_width(s), solver_field_height(s));
421 }
422
solver_set_field(solver_t * s,starxy_t * field)423 void solver_set_field(solver_t* s, starxy_t* field) {
424 if (s->fieldxy)
425 starxy_free(s->fieldxy);
426 s->fieldxy = field;
427 // Preprocessing happens in "solver_preprocess_field()".
428 }
429
solver_set_field_bounds(solver_t * s,double xlo,double xhi,double ylo,double yhi)430 void solver_set_field_bounds(solver_t* s, double xlo, double xhi, double ylo, double yhi) {
431 s->field_minx = xlo;
432 s->field_maxx = xhi;
433 s->field_miny = ylo;
434 s->field_maxy = yhi;
435 set_diag(s);
436 }
437
solver_cleanup_field(solver_t * solver)438 void solver_cleanup_field(solver_t* solver) {
439 solver_reset_best_match(solver);
440 solver_free_field(solver);
441 solver->fieldxy = NULL;
442 solver_reset_counters(solver);
443 }
444
solver_verify_sip_wcs(solver_t * solver,sip_t * sip)445 void solver_verify_sip_wcs(solver_t* solver, sip_t* sip) { //, MatchObj* pmo) {
446 int i, nindexes;
447 MatchObj mo;
448 MatchObj* pmo;
449 anbool olddqb;
450
451 pmo = &mo;
452
453 if (!solver->vf)
454 solver_preprocess_field(solver);
455
456 // fabricate a match and inject it into the solver.
457 set_matchobj_template(solver, pmo);
458 memcpy(&(mo.wcstan), &(sip->wcstan), sizeof(tan_t));
459 mo.wcs_valid = TRUE;
460 mo.scale = sip_pixel_scale(sip);
461 set_center_and_radius(solver, pmo, NULL, sip);
462 olddqb = solver->distance_from_quad_bonus;
463 solver->distance_from_quad_bonus = FALSE;
464
465 nindexes = pl_size(solver->indexes);
466 for (i=0; i<nindexes; i++) {
467 index_t* index = pl_get(solver->indexes, i);
468 set_index(solver, index);
469 solver_inject_match(solver, pmo, sip);
470 }
471
472 // revert
473 solver->distance_from_quad_bonus = olddqb;
474 }
475
solver_add_index(solver_t * solver,index_t * index)476 void solver_add_index(solver_t* solver, index_t* index) {
477 pl_append(solver->indexes, index);
478 }
479
solver_n_indices(const solver_t * solver)480 int solver_n_indices(const solver_t* solver) {
481 return pl_size(solver->indexes);
482 }
483
solver_get_index(const solver_t * solver,int i)484 index_t* solver_get_index(const solver_t* solver, int i) {
485 return pl_get(solver->indexes, i);
486 }
487
solver_reset_best_match(solver_t * sp)488 void solver_reset_best_match(solver_t* sp) {
489 // we don't really care about very bad best matches...
490 sp->best_logodds = 0;
491 memset(&(sp->best_match), 0, sizeof(MatchObj));
492 sp->best_index = NULL;
493 sp->best_match_solves = FALSE;
494 sp->have_best_match = FALSE;
495 }
496
solver_compute_quad_range(const solver_t * sp,const index_t * index,double * minAB,double * maxAB)497 void solver_compute_quad_range(const solver_t* sp, const index_t* index,
498 double* minAB, double* maxAB) {
499 double scalefudge; // in pixels
500
501 // compute fudge factor for quad scale: what are the extreme
502 // ranges of quad scales that should be accepted, given the
503 // code tolerance?
504 // -what is the maximum number of pixels a C or D star can move
505 // to singlehandedly exceed the code tolerance?
506 // -largest quad
507 // -smallest arcsec-per-pixel scale
508
509 // -index_scale_upper * 1/sqrt(2) is the side length of
510 // the unit-square of code space, in arcseconds.
511 // -that times the code tolerance is how far a C/D star
512 // can move before exceeding the code tolerance, in arcsec.
513 // -that divided by the smallest arcsec-per-pixel scale
514 // gives the largest motion in pixels.
515
516 //logverb("Index scale %f, %f\n",
517 //index->index_scale_upper, index->index_scale_lower);
518
519 scalefudge = index->index_scale_upper * M_SQRT1_2 *
520 sp->codetol / sp->funits_upper;
521
522 if (sp->funits_upper != 0.0) {
523 *minAB = index->index_scale_lower / sp->funits_upper;
524 *minAB -= scalefudge;
525 }
526 if (sp->funits_lower != 0.0) {
527 *maxAB = index->index_scale_upper / sp->funits_lower;
528 *maxAB += scalefudge;
529 }
530 }
531
532 static void try_all_codes(const pquad* pq,
533 const int* fieldstars, int dimquad,
534 solver_t* solver, double tol2);
535
536 static void try_all_codes_2(const int* fieldstars, int dimquad,
537 const double* code, solver_t* solver,
538 anbool current_parity, double tol2);
539
540 static void try_permutations(const int* origstars, int dimquad,
541 const double* origcode,
542 solver_t* solver, anbool current_parity,
543 double tol2,
544 int* stars, double* code,
545 int slot, anbool* placed,
546 kdtree_qres_t** presult);
547
548 static void resolve_matches(kdtree_qres_t* krez, const double *field,
549 const int* fstars, int dimquads,
550 solver_t* solver, anbool current_parity);
551
552 static int solver_handle_hit(solver_t* sp, MatchObj* mo, sip_t* sip, anbool fake_match);
553
check_scale(pquad * pq,solver_t * s)554 static void check_scale(pquad* pq, solver_t* s) {
555 double dx, dy;
556 dx = field_getx(s, pq->fieldB) - field_getx(s, pq->fieldA);
557 dy = field_gety(s, pq->fieldB) - field_gety(s, pq->fieldA);
558 pq->scale = dx*dx + dy*dy;
559 if ((pq->scale < s->minminAB2) ||
560 (pq->scale > s->maxmaxAB2)) {
561 pq->scale_ok = FALSE;
562 return;
563 }
564 pq->costheta = (dy + dx) / pq->scale;
565 pq->sintheta = (dy - dx) / pq->scale;
566 pq->rel_field_noise2 = (s->verify_pix * s->verify_pix) / pq->scale;
567 pq->scale_ok = TRUE;
568 }
569
check_inbox(pquad * pq,int start,solver_t * solver)570 static void check_inbox(pquad* pq, int start, solver_t* solver) {
571 int i;
572 double Ax, Ay;
573 field_getxy(solver, pq->fieldA, &Ax, &Ay);
574 // check which C, D points are inside the circle.
575 for (i = start; i < pq->ninbox; i++) {
576 double r;
577 double Cx, Cy, xxtmp;
578 double tol = solver->codetol;
579 if (!pq->inbox[i])
580 continue;
581 field_getxy(solver, i, &Cx, &Cy);
582 Cx -= Ax;
583 Cy -= Ay;
584 xxtmp = Cx;
585 Cx = Cx * pq->costheta + Cy * pq->sintheta;
586 Cy = -xxtmp * pq->sintheta + Cy * pq->costheta;
587
588 // make sure it's in the circle centered at (0.5, 0.5)
589 // with radius 1/sqrt(2) (plus codetol for fudge):
590 // (x-1/2)^2 + (y-1/2)^2 <= (r + codetol)^2
591 // x^2-x+1/4 + y^2-y+1/4 <= (1/sqrt(2) + codetol)^2
592 // x^2-x + y^2-y + 1/2 <= 1/2 + sqrt(2)*codetol + codetol^2
593 // x^2-x + y^2-y <= sqrt(2)*codetol + codetol^2
594 r = (Cx * Cx - Cx) + (Cy * Cy - Cy);
595 if (r > (tol * (M_SQRT2 + tol))) {
596 pq->inbox[i] = FALSE;
597 continue;
598 }
599 setx(pq->xy, i, Cx);
600 sety(pq->xy, i, Cy);
601 }
602 }
603
604 #if defined DEBUGSOLVER
print_inbox(pquad * pq)605 static void print_inbox(pquad* pq) {
606 int i;
607 debug("[ ");
608 for (i = 0; i < pq->ninbox; i++) {
609 if (pq->inbox[i])
610 debug("%i ", i);
611 }
612 debug("] (n %i)\n", pq->ninbox);
613 }
614 #else
print_inbox(pquad * pq)615 static void print_inbox(pquad* pq) {}
616 #endif
617
618
solver_reset_field_size(solver_t * s)619 void solver_reset_field_size(solver_t* s) {
620 s->field_minx = s->field_maxx = s->field_miny = s->field_maxy = 0;
621 s->field_diag = 0.0;
622 }
623
find_field_boundaries(solver_t * solver)624 static void find_field_boundaries(solver_t* solver) {
625 // If the bounds haven't been set, use the bounding box.
626 if ((solver->field_minx == solver->field_maxx) ||
627 (solver->field_miny == solver->field_maxy)) {
628 int i;
629 solver->field_minx = solver->field_miny = HUGE_VAL;
630 solver->field_maxx = solver->field_maxy = -HUGE_VAL;
631 for (i = 0; i < starxy_n(solver->fieldxy); i++) {
632 solver->field_minx = MIN(solver->field_minx, field_getx(solver, i));
633 solver->field_maxx = MAX(solver->field_maxx, field_getx(solver, i));
634 solver->field_miny = MIN(solver->field_miny, field_gety(solver, i));
635 solver->field_maxy = MAX(solver->field_maxy, field_gety(solver, i));
636 }
637 }
638 set_diag(solver);
639 }
640
solver_preprocess_field(solver_t * solver)641 void solver_preprocess_field(solver_t* solver) {
642 find_field_boundaries(solver);
643 // precompute a kdtree over the field
644 solver->vf = verify_field_preprocess(solver->fieldxy);
645
646 solver->vf->do_uniformize = solver->verify_uniformize;
647 solver->vf->do_dedup = solver->verify_dedup;
648 }
649
solver_free_field(solver_t * solver)650 void solver_free_field(solver_t* solver) {
651 //# Modified by Robert Lancaster for the StellarSolver Internal Library
652 //if (solver->fieldxy)
653 // starxy_free(solver->fieldxy);
654 //solver->fieldxy = NULL;
655 if (solver->vf)
656 verify_field_free(solver->vf);
657 solver->vf = NULL;
658 }
659
solver_get_field(solver_t * solver)660 starxy_t* solver_get_field(solver_t* solver) {
661 return solver->fieldxy;
662 }
663
get_tolerance(solver_t * solver)664 static double get_tolerance(solver_t* solver) {
665 return square(solver->codetol);
666 /*
667 double maxtol2 = square(solver->codetol);
668 double tol2;
669 tol2 = 49.0 * (solver->rel_field_noise2 + solver->rel_index_noise2);
670 //printf("code tolerance %g.\n", sqrt(tol2));
671 if (tol2 > maxtol2)
672 tol2 = maxtol2;
673 return tol2;
674 */
675 }
676
677 /*
678 A somewhat tricky recursive function: stars A and B have already been
679 chosen, so the code coordinate system has been fixed, and we've
680 already determined which other stars will create valid codes (ie, are
681 in the "box"). Now we want to build features using all sets of valid
682 stars (without permutations).
683
684 pq - data associated with the AB pair.
685 field - the array of field star numbers
686 fieldoffset - offset into the field array where we should add the first star
687 n_to_add - number of stars to add
688 adding - the star we're currently adding; in [0, n_to_add).
689 fieldtop - the maximum field star number to build quads out of.
690 dimquad, solver, tol2 - passed to try_all_codes.
691 */
add_stars(const pquad * pq,int * field,int fieldoffset,int n_to_add,int adding,int fieldtop,int dimquad,solver_t * solver,double tol2)692 static void add_stars(const pquad* pq, int* field, int fieldoffset,
693 int n_to_add, int adding, int fieldtop,
694 int dimquad,
695 solver_t* solver, double tol2) {
696 int bottom;
697 int* f = field + fieldoffset;
698 // When we're adding the first star, we start from index zero.
699 // When we're adding subsequent stars, we start from the previous value
700 // plus one, to avoid adding permutations.
701 bottom = (adding ? f[adding-1] + 1 : 0);
702
703 // It looks funny that we're using f[adding] as a loop variable, but
704 // it's required because try_all_codes needs to know which field stars
705 // were used to create the quad (which are stored in the "f" array)
706 for (f[adding]=bottom; f[adding]<fieldtop; f[adding]++) {
707 if (!pq->inbox[f[adding]])
708 continue;
709 if (unlikely(solver->quit_now))
710 return;
711
712 // If we've hit the end of the recursion (we're adding the last star),
713 // call try_all_codes to try the quad we've built.
714 if (adding == n_to_add-1) {
715 // (when not testing, TRY_ALL_CODES is just try_all_codes.)
716 TRY_ALL_CODES(pq, field, dimquad, solver, tol2);
717 } else {
718 // Else recurse.
719 add_stars(pq, field, fieldoffset, n_to_add, adding+1,
720 fieldtop, dimquad, solver, tol2);
721 }
722 }
723 }
724
725
726 // The real deal
solver_run(solver_t * solver)727 void solver_run(solver_t* solver) {
728 int numxy, newpoint;
729 double usertime, systime;
730 // first timer callback is called after 1 second
731 time_t next_timer_callback_time = time(NULL) + 1;
732 pquad* pquads;
733 size_t i, num_indexes;
734 double tol2;
735 int field[DQMAX];
736
737 get_resource_stats(&usertime, &systime, NULL);
738
739 if (!solver->vf)
740 solver_preprocess_field(solver);
741
742 memset(field, 0, sizeof(field));
743
744 solver->starttime = usertime + systime;
745
746 numxy = starxy_n(solver->fieldxy);
747 if (solver->endobj && (numxy > solver->endobj))
748 numxy = solver->endobj;
749 if (solver->startobj >= numxy)
750 return;
751 if (numxy >= 1000) {
752 logverb("Limiting search to first 1000 objects\n");
753 numxy = 1000;
754 }
755
756 if (solver->set_crpix && solver->set_crpix_center) {
757 solver->crpix[0] = wcs_pixel_center_for_size(solver_field_width(solver));
758 solver->crpix[1] = wcs_pixel_center_for_size(solver_field_height(solver));
759 logverb("Setting CRPIX to center (%.1f, %.1f) based on image size %i x %i\n",
760 solver->crpix[0], solver->crpix[1],
761 (int)solver_field_width(solver), (int)solver_field_height(solver));
762 }
763
764 num_indexes = pl_size(solver->indexes);
765 {
766 #ifndef _MSC_VER //# Modified by Robert Lancaster for the StellarSolver Internal Library
767 double minAB2s[num_indexes];
768 double maxAB2s[num_indexes];
769 #else
770 double* minAB2s = (double*)malloc(sizeof(double)*num_indexes);
771 double* maxAB2s = (double*)malloc(sizeof(double)*num_indexes);
772 #endif
773 solver->minminAB2 = HUGE_VAL;
774 solver->maxmaxAB2 = -HUGE_VAL;
775 for (i = 0; i < num_indexes; i++) {
776 double minAB=0, maxAB=0;
777 index_t* index = pl_get(solver->indexes, i);
778 // The limits on the size of quads that we try to match, in pixels.
779 // Derived from index_scale_* and funits_*.
780 solver_compute_quad_range(solver, index, &minAB, &maxAB);
781 //logverb("Index \"%s\" quad range %f to %f\n", index->indexname,
782 //minAB, maxAB);
783 minAB2s[i] = square(minAB);
784 maxAB2s[i] = square(maxAB);
785 solver->minminAB2 = MIN(solver->minminAB2, minAB2s[i]);
786 solver->maxmaxAB2 = MAX(solver->maxmaxAB2, maxAB2s[i]);
787
788 if (index->cx_less_than_dx) {
789 solver->cxdx_margin = 1.5 * solver->codetol;
790 // FIXME die horribly if the indexes have differing cx_less_than_dx
791 }
792 }
793 solver->minminAB2 = MAX(solver->minminAB2, square(solver->quadsize_min));
794 if (solver->quadsize_max != 0.0)
795 solver->maxmaxAB2 = MIN(solver->maxmaxAB2, square(solver->quadsize_max));
796 logverb("Quad scale range: [%g, %g] pixels\n", sqrt(solver->minminAB2), sqrt(solver->maxmaxAB2));
797
798 // quick-n-dirty scale estimate using stars A,B.
799 solver->abscale_high = square(arcsec2rad(solver->funits_upper) * (1.0 + solver->codetol));
800 solver->abscale_low = square(arcsec2rad(solver->funits_lower) * (1.0 - solver->codetol));
801
802 /** Ugh, I want to avoid doing distsq2rad when checking scale,
803 but that means correcting for the difference between the
804 distance along the curve of the sphere vs the chord distance.
805 This affects the lower bound for the largest quads in a messy way...
806 This below isn't right.
807 solver->abscale_high = square(arcsec2rad(solver->funits_upper) * (1.0 + solver->codetol));
808 solver->abscale_low = arcsec2rad(solver->funits_lower) * (1.0 - solver->codetol) *
809 MIN(M_PI, arcsec2rad(field_diag * solver->funits_upper)) ...
810 */
811
812 pquads = calloc(numxy * numxy, sizeof(pquad));
813
814 /* We maintain an array of "potential quads" (pquad) structs, where
815 * each struct corresponds to one choice of stars A and B; the struct
816 * at index (B * numxy + A) holds information about quads that could be
817 * created using stars A,B.
818 *
819 * (We only use the above-diagonal elements of this 2D array because
820 * A<B.)
821 *
822 * For each AB pair, we cache the scale and the rotation parameters,
823 * and we keep an array "inbox" of length "numxy" of booleans, one for
824 * each star, which say whether that star is eligible to be star C or D
825 * of a quad with AB at the corners. (Obviously A and B aren't
826 * eligible).
827 *
828 * The "ninbox" parameter is somewhat misnamed - it says that "inbox"
829 * elements in the range [0, ninbox) have been initialized.
830 */
831
832 /* (See explanatory paragraph below) If "solver->startobj" isn't zero,
833 * then we need to initialize the triangle of "pquads" up to
834 * A=startobj-2, B=startobj-1. */
835 if (solver->startobj) {
836 debug("startobj > 0; priming pquad arrays.\n");
837 for (field[B] = 0; field[B] < solver->startobj; field[B]++) {
838 for (field[A] = 0; field[A] < field[B]; field[A]++) {
839 pquad* pq = pquads + field[B] * numxy + field[A];
840 pq->fieldA = field[A];
841 pq->fieldB = field[B];
842 debug("trying A=%i, B=%i\n", field[A], field[B]);
843 check_scale(pq, solver);
844 if (!pq->scale_ok) {
845 debug(" bad scale for A=%i, B=%i\n", field[A], field[B]);
846 continue;
847 }
848 pq->xy = malloc(numxy * 2 * sizeof(double));
849 pq->inbox = malloc(numxy * sizeof(anbool));
850 memset(pq->inbox, TRUE, solver->startobj);
851 pq->ninbox = solver->startobj;
852 pq->inbox[field[A]] = FALSE;
853 pq->inbox[field[B]] = FALSE;
854 check_inbox(pq, 0, solver);
855 debug(" inbox(A=%i, B=%i): ", field[A], field[B]);
856 print_inbox(pq);
857 }
858 }
859 }
860
861 /* Each time through the "for" loop below, we consider a new star
862 * ("newpoint"). First, we try building all quads that have the new
863 * star on the diagonal (star B). Then, we try building all quads that
864 * have the star not on the diagonal (star D).
865 *
866 * For each AB pair, we have a "potential_quad" or "pquad" struct.
867 * This caches the computation we need to do: deciding whether the
868 * scale is acceptable, computing the transformation to code
869 * coordinates, and deciding which C,D stars are in the circle.
870 */
871 for (newpoint = solver->startobj; newpoint < numxy; newpoint++) {
872
873 debug("Trying newpoint=%i\n", newpoint);
874
875 // Give our caller a chance to cancel us midway. The callback
876 // returns how long to wait before calling again.
877
878 if (solver->timer_callback) {
879 time_t delay;
880 time_t now = time(NULL);
881 if (now > next_timer_callback_time) {
882 update_timeused(solver);
883 delay = solver->timer_callback(solver->userdata);
884 if (delay == 0) // Canceled
885 break;
886 next_timer_callback_time = now + delay;
887 }
888 }
889
890 solver->last_examined_object = newpoint;
891 // quads with the new star on the diagonal:
892 field[B] = newpoint;
893 debug("Trying quads with B=%i\n", newpoint);
894
895 // first do an index-independent scale check...
896 for (field[A] = 0; field[A] < newpoint; field[A]++) {
897 // initialize the "pquad" struct for this AB combo.
898 pquad* pq = pquads + field[B] * numxy + field[A];
899 pq->fieldA = field[A];
900 pq->fieldB = field[B];
901 debug(" trying A=%i, B=%i\n", field[A], field[B]);
902 check_scale(pq, solver);
903 if (!pq->scale_ok) {
904 debug(" bad scale for A=%i, B=%i\n", field[A], field[B]);
905 continue;
906 }
907 // initialize the "inbox" array:
908 pq->inbox = malloc(numxy * sizeof(anbool));
909 pq->xy = malloc(numxy * 2 * sizeof(double));
910 // -try all stars up to "newpoint"...
911 assert(sizeof(anbool) == 1);
912 memset(pq->inbox, TRUE, newpoint + 1);
913 pq->ninbox = newpoint + 1;
914 // -except A and B.
915 pq->inbox[field[A]] = FALSE;
916 pq->inbox[field[B]] = FALSE;
917 check_inbox(pq, 0, solver);
918 debug(" inbox(A=%i, B=%i): ", field[A], field[B]);
919 print_inbox(pq);
920 }
921
922 // Now iterate through the different indices
923 for (i = 0; i < num_indexes; i++) {
924 index_t* index = pl_get(solver->indexes, i);
925 int dimquads;
926 set_index(solver, index);
927 dimquads = index_dimquads(index);
928 for (field[A] = 0; field[A] < newpoint; field[A]++) {
929 // initialize the "pquad" struct for this AB combo.
930 pquad* pq = pquads + field[B] * numxy + field[A];
931 if (!pq->scale_ok)
932 continue;
933 if ((pq->scale < minAB2s[i]) ||
934 (pq->scale > maxAB2s[i]))
935 continue;
936 // set code tolerance for this index and AB pair...
937 solver->rel_field_noise2 = pq->rel_field_noise2;
938 tol2 = get_tolerance(solver);
939 // Now look at all sets of (C, D, ...) stars (subject to field[C] < field[D] < ...)
940 // ("dimquads - 2" because we've set stars A and B at this point)
941 add_stars(pq, field, C, dimquads-2, 0, newpoint, dimquads, solver, tol2);
942 if (solver->quit_now)
943 goto quitnow;
944 }
945 }
946
947 if (solver->quit_now)
948 goto quitnow;
949
950 // Now try building quads with the new star not on the diagonal:
951 field[C] = newpoint;
952 // (in this loop field[C] > field[D])
953 debug("Trying quads with C=%i\n", newpoint);
954 for (field[A] = 0; field[A] < newpoint; field[A]++) {
955 for (field[B] = field[A] + 1; field[B] < newpoint; field[B]++) {
956 // grab the "pquad" for this AB combo
957 pquad* pq = pquads + field[B] * numxy + field[A];
958 if (!pq->scale_ok) {
959 debug(" bad scale for A=%i, B=%i\n", field[A], field[B]);
960 continue;
961 }
962 // test if this C is in the box:
963 pq->inbox[field[C]] = TRUE;
964 pq->ninbox = field[C] + 1;
965 check_inbox(pq, field[C], solver);
966 if (!pq->inbox[field[C]]) {
967 debug(" C is not in the box for A=%i, B=%i\n", field[A], field[B]);
968 continue;
969 }
970 debug(" C is in the box for A=%i, B=%i\n", field[A], field[B]);
971 debug(" box now:");
972 print_inbox(pq);
973 debug("\n");
974
975 solver->rel_field_noise2 = pq->rel_field_noise2;
976
977 for (i = 0; i < pl_size(solver->indexes); i++) {
978 int dimquads;
979 index_t* index = pl_get(solver->indexes, i);
980 if ((pq->scale < minAB2s[i]) ||
981 (pq->scale > maxAB2s[i]))
982 continue;
983 set_index(solver, index);
984 dimquads = index_dimquads(index);
985
986 tol2 = get_tolerance(solver);
987
988 if (dimquads > 3) {
989 // ("dimquads - 3" because we've set stars A, B, and C at this point)
990 add_stars(pq, field, D, dimquads-3, 0, newpoint, dimquads, solver, tol2);
991 } else {
992 TRY_ALL_CODES(pq, field, dimquads, solver, tol2);
993 }
994 if (solver->quit_now)
995 goto quitnow;
996 }
997 }
998 }
999 logverb("object %u of %u: %i quads tried, %i matched.\n",
1000 newpoint + 1, numxy, solver->numtries, solver->nummatches);
1001
1002 if ((solver->maxquads && (solver->numtries >= solver->maxquads))
1003 || (solver->maxmatches && (solver->nummatches >= solver->maxmatches))
1004 || solver->quit_now)
1005 break;
1006 }
1007
1008 quitnow:
1009 for (i = 0; i < (numxy*numxy); i++) {
1010 pquad* pq = pquads + i;
1011 free(pq->inbox);
1012 free(pq->xy);
1013 }
1014 free(pquads);
1015
1016 #ifdef _MSC_VER //# Modified by Robert Lancaster for the StellarSolver Internal Library
1017 free(minAB2s);
1018 free(maxAB2s);
1019 #endif
1020 }
1021 }
1022
1023 /**
1024 All the stars in this quad have been chosen. Figure out which
1025 permutations of stars CDE are valid and search for matches.
1026 */
try_all_codes(const pquad * pq,const int * fieldstars,int dimquad,solver_t * solver,double tol2)1027 static void try_all_codes(const pquad* pq,
1028 const int* fieldstars, int dimquad,
1029 solver_t* solver, double tol2) {
1030 int dimcode = (dimquad - 2) * 2;
1031 double code[DCMAX];
1032 double flipcode[DCMAX];
1033 int i;
1034
1035 solver->numtries++;
1036
1037 debug(" trying quad [");
1038 for (i=0; i<dimquad; i++) {
1039 debug("%s%i", (i?" ":""), fieldstars[i]);
1040 }
1041 debug("]\n");
1042
1043 for (i=0; i<dimquad-NBACK; i++) {
1044 code[2*i ] = getx(pq->xy, fieldstars[NBACK+i]);
1045 code[2*i+1] = gety(pq->xy, fieldstars[NBACK+i]);
1046 }
1047
1048 if (solver->parity == PARITY_NORMAL ||
1049 solver->parity == PARITY_BOTH) {
1050
1051 debug(" trying normal parity: code=[");
1052 for (i=0; i<dimcode; i++)
1053 debug("%s%g", (i?", ":""), code[i]);
1054 debug("].\n");
1055
1056 try_all_codes_2(fieldstars, dimquad, code, solver, FALSE, tol2);
1057 }
1058 if (solver->parity == PARITY_FLIP ||
1059 solver->parity == PARITY_BOTH) {
1060
1061 quad_flip_parity(code, flipcode, dimcode);
1062
1063 debug(" trying reverse parity: code=[");
1064 for (i=0; i<dimcode; i++)
1065 debug("%s%g", (i?", ":""), flipcode[i]);
1066 debug("].\n");
1067
1068 try_all_codes_2(fieldstars, dimquad, flipcode, solver, TRUE, tol2);
1069 }
1070 }
1071
1072 /**
1073 This function tries the quad with the "backbone" stars A and B in
1074 normal and flipped configurations.
1075 */
try_all_codes_2(const int * fieldstars,int dimquad,const double * code,solver_t * solver,anbool current_parity,double tol2)1076 static void try_all_codes_2(const int* fieldstars, int dimquad,
1077 const double* code, solver_t* solver,
1078 anbool current_parity, double tol2) {
1079 int i;
1080 kdtree_qres_t* result = NULL;
1081 int dimcode = (dimquad - 2) * 2;
1082 int stars[DQMAX];
1083 double flipcode[DCMAX];
1084
1085 // We actually only use elements up to dimquads-2.
1086 anbool placed[DQMAX];
1087
1088 // Un-flipped:
1089 stars[0] = fieldstars[0];
1090 stars[1] = fieldstars[1];
1091
1092 for (i=0; i<DQMAX; i++)
1093 placed[i] = FALSE;
1094
1095 try_permutations(fieldstars, dimquad, code, solver, current_parity,
1096 tol2, stars, NULL, 0, placed, &result);
1097 if (unlikely(solver->quit_now))
1098 goto bailout;
1099
1100 // Flipped:
1101 stars[0] = fieldstars[1];
1102 stars[1] = fieldstars[0];
1103
1104 for (i=0; i<dimcode; i++)
1105 flipcode[i] = 1.0 - code[i];
1106
1107 for (i=0; i<DQMAX; i++)
1108 placed[i] = FALSE;
1109
1110 try_permutations(fieldstars, dimquad, flipcode, solver, current_parity,
1111 tol2, stars, NULL, 0, placed, &result);
1112
1113 bailout:
1114 kdtree_free_query(result);
1115 }
1116
1117 /**
1118 This functions tries different permutations of the non-backbone
1119 stars C [, D [,E ] ]
1120 */
try_permutations(const int * origstars,int dimquad,const double * origcode,solver_t * solver,anbool current_parity,double tol2,int * stars,double * code,int slot,anbool * placed,kdtree_qres_t ** presult)1121 static void try_permutations(const int* origstars, int dimquad,
1122 const double* origcode,
1123 solver_t* solver, anbool current_parity,
1124 double tol2,
1125 int* stars, double* code,
1126 int slot, anbool* placed,
1127 kdtree_qres_t** presult) {
1128 int i;
1129 int options = KD_OPTIONS_SMALL_RADIUS | KD_OPTIONS_COMPUTE_DISTS |
1130 KD_OPTIONS_NO_RESIZE_RESULTS | KD_OPTIONS_USE_SPLIT;
1131 double mycode[DCMAX];
1132 int Nstars = dimquad - NBACK;
1133 int lastslot = dimquad - NBACK - 1;
1134 /*
1135 This is a recursive function that tries all combinations of the
1136 "internal" stars (ie, not stars A,B that form the "backbone" of
1137 the quad).
1138
1139 We fill the "stars" array with the star IDs (from "origstars") of
1140 the stars that form the quad, while simultaneously filling the
1141 "code" array with the corresponding code coordinates (from
1142 "origcode").
1143
1144 For example, if "dimquad" is 5, and "origstars" contains
1145 A,B,C,D,E, we want to call "resolve_matches" with the following
1146 combinations in "stars":
1147
1148 AB CDE
1149 AB CED
1150 AB DCE
1151 AB DEC
1152 AB ECD
1153 AB EDC
1154
1155 This call will try to put each star in "slot" in turn, then for
1156 each one recurse to "slot" in the rest of the stars.
1157
1158 Note that we are filling stars[2], stars[3], etc; the first two
1159 elements are already filled by stars A and B.
1160 */
1161
1162 if (code == NULL)
1163 code = mycode;
1164
1165 // We try putting each star that hasn't already been placed in
1166 // this "slot".
1167 for (i=0; i<Nstars; i++) {
1168 if (placed[i])
1169 continue;
1170
1171 // Check cx <= dx, if we're a "dx".
1172 if (slot > 0 && solver->index->cx_less_than_dx) {
1173 if (code[2*(slot - 1) +0] > origcode[2*i +0] + solver->cxdx_margin) {
1174 debug("cx <= dx check failed: %g > %g + %g\n",
1175 code[2*(slot - 1) +0], origcode[2*i +0],
1176 solver->cxdx_margin);
1177 solver->num_cxdx_skipped++;
1178 continue;
1179 }
1180 }
1181
1182 // Slot in this star...
1183 stars[slot + NBACK] = origstars[i + NBACK];
1184 code[2*slot +0] = origcode[2*i +0];
1185 code[2*slot +1] = origcode[2*i +1];
1186
1187 // Check meanx <= 1/2.
1188 if (solver->index->cx_less_than_dx &&
1189 solver->index->meanx_less_than_half) {
1190 // Check the "cx + dx <= 1" condition (for quads); in general,
1191 // combined with the "cx <= dx" condition, this means that the
1192 // mean(x) <= 1/2.
1193 int j;
1194 double meanx = 0;
1195 for (j=0; j<=slot; j++)
1196 meanx += code[2*j];
1197 meanx /= (slot+1);
1198 if (meanx > 0.5 + solver->cxdx_margin) {
1199 debug("meanx <= 0.5 check failed: %g > 0.5 + %g\n",
1200 meanx, solver->cxdx_margin);
1201 solver->num_meanx_skipped++;
1202 continue;
1203 }
1204 }
1205
1206 // If we have more slots to fill...
1207 if (slot < lastslot) {
1208 placed[i] = TRUE;
1209 try_permutations(origstars, dimquad, origcode, solver,
1210 current_parity, tol2, stars, code,
1211 slot+1, placed, presult);
1212 placed[i] = FALSE;
1213
1214 } else {
1215 #if defined(TESTING_TRYPERMUTATIONS)
1216 TEST_TRY_PERMUTATIONS(stars, code, dimquad, solver);
1217 continue;
1218 #endif
1219
1220 // Search with the code we've built.
1221 *presult = kdtree_rangesearch_options_reuse
1222 (solver->index->codekd->tree, *presult, code, tol2, options);
1223 //debug(" trying ABCD = [%i %i %i %i]: %i results.\n",
1224 //fstars[A], fstars[B], fstars[C], fstars[D], result->nres);
1225
1226 if ((*presult)->nres) {
1227 double pixvals[DQMAX*2];
1228 int j;
1229 for (j=0; j<dimquad; j++) {
1230 setx(pixvals, j, field_getx(solver, stars[j]));
1231 sety(pixvals, j, field_gety(solver, stars[j]));
1232 }
1233 resolve_matches(*presult, pixvals, stars, dimquad, solver,
1234 current_parity);
1235 }
1236 if (unlikely(solver->quit_now))
1237 return;
1238 }
1239 }
1240 }
1241
1242 // "field" contains the xy pixel coordinates of stars A,B,C,D.
resolve_matches(kdtree_qres_t * krez,const double * field,const int * fieldstars,int dimquads,solver_t * solver,anbool current_parity)1243 static void resolve_matches(kdtree_qres_t* krez, const double *field,
1244 const int* fieldstars, int dimquads,
1245 solver_t* solver, anbool current_parity) {
1246 int jj, thisquadno;
1247 MatchObj mo;
1248
1249 #ifndef _MSC_VER //# Modified by Robert Lancaster for the StellarSolver Internal Library
1250 unsigned int star[dimquads];
1251 #else
1252 unsigned int* star = (unsigned int*)malloc(sizeof(unsigned int)*dimquads);
1253 #endif
1254
1255 assert(krez);
1256
1257 for (jj = 0; jj < krez->nres; jj++) {
1258
1259 #ifndef _MSC_VER //# Modified by Robert Lancaster for the StellarSolver Internal Library
1260 double starxyz[dimquads*3];
1261 #else
1262 double* starxyz = (double*)malloc(sizeof(double)*dimquads*3);
1263 #endif
1264
1265 double scale;
1266 double arcsecperpix;
1267 tan_t wcs;
1268 int i;
1269 anbool outofbounds = FALSE;
1270 double abscale;
1271
1272 solver->nummatches++;
1273 thisquadno = krez->inds[jj];
1274 quadfile_get_stars(solver->index->quads, thisquadno, star);
1275 for (i=0; i<dimquads; i++) {
1276 startree_get(solver->index->starkd, star[i], starxyz + 3*i);
1277 if (solver->use_radec)
1278 if (distsq(starxyz + 3*i, solver->centerxyz, 3) > solver->r2) {
1279 outofbounds = TRUE;
1280 break;
1281 }
1282 }
1283 if (outofbounds) {
1284 debug("Quad match is out of bounds.\n");
1285 solver->num_radec_skipped++;
1286 continue;
1287 }
1288
1289 debug(" stars [");
1290 for (i=0; i<dimquads; i++)
1291 debug("%s%i", (i?" ":""), star[i]);
1292 debug("]\n");
1293
1294 // Quick-n-dirty scale estimate based on two stars.
1295 //abscale = distsq(starxyz, starxyz+3, 3) / distsq(field, field+2, 2);
1296 // in (rad per pix)**2
1297 abscale = square(distsq2rad(distsq(starxyz, starxyz+3, 3))) /
1298 distsq(field, field+2, 2);
1299 if (abscale > solver->abscale_high ||
1300 abscale < solver->abscale_low) {
1301 solver->num_abscale_skipped++;
1302 continue;
1303 }
1304
1305 // compute TAN projection from the matching quad alone.
1306 if (fit_tan_wcs(starxyz, field, dimquads, &wcs, &scale)) {
1307 // bad quad.
1308 logverb("bad quad at %s:%i\n", __FILE__, __LINE__);
1309 continue;
1310 }
1311 arcsecperpix = scale * 3600.0;
1312
1313 // FIXME - should there be scale fudge here?
1314 if (arcsecperpix > solver->funits_upper ||
1315 arcsecperpix < solver->funits_lower) {
1316 debug(" bad scale (%g arcsec/pix, range %g %g)\n",
1317 arcsecperpix, solver->funits_lower, solver->funits_upper);
1318 continue;
1319 }
1320 solver->numscaleok++;
1321
1322 set_matchobj_template(solver, &mo);
1323 memcpy(&(mo.wcstan), &wcs, sizeof(tan_t));
1324 mo.wcs_valid = TRUE;
1325 mo.code_err = krez->sdists[jj];
1326 mo.scale = arcsecperpix;
1327 mo.parity = current_parity;
1328 mo.quads_tried = solver->numtries;
1329 mo.quads_matched = solver->nummatches;
1330 mo.quads_scaleok = solver->numscaleok;
1331 mo.quad_npeers = krez->nres;
1332 mo.timeused = solver->timeused;
1333 mo.quadno = thisquadno;
1334 mo.dimquads = dimquads;
1335 for (i=0; i<dimquads; i++) {
1336 mo.star[i] = star[i];
1337 mo.field[i] = fieldstars[i];
1338 mo.ids[i] = 0;
1339 }
1340
1341 memcpy(mo.quadpix, field, 2 * dimquads * sizeof(double));
1342 memcpy(mo.quadxyz, starxyz, 3 * dimquads * sizeof(double));
1343
1344 set_center_and_radius(solver, &mo, &(mo.wcstan), NULL);
1345
1346 if (solver_handle_hit(solver, &mo, NULL, FALSE))
1347 solver->quit_now = TRUE;
1348
1349 if (unlikely(solver->quit_now))
1350 {
1351 #ifdef _MSC_VER //# Modified by Robert Lancaster for the StellarSolver Internal Library
1352 free(starxyz);
1353 free(star);
1354 #endif
1355 return;
1356 }
1357
1358 #ifdef _MSC_VER //# Modified by Robert Lancaster for the StellarSolver Internal Library
1359 free(starxyz);
1360 #endif
1361 }
1362 #ifdef _MSC_VER //# Modified by Robert Lancaster for the StellarSolver Internal Library
1363 free(star);
1364 #endif
1365 }
1366
solver_inject_match(solver_t * solver,MatchObj * mo,sip_t * sip)1367 void solver_inject_match(solver_t* solver, MatchObj* mo, sip_t* sip) {
1368 solver_handle_hit(solver, mo, sip, TRUE);
1369 }
1370
solver_handle_hit(solver_t * sp,MatchObj * mo,sip_t * sip,anbool fake_match)1371 static int solver_handle_hit(solver_t* sp, MatchObj* mo, sip_t* sip,
1372 anbool fake_match) {
1373 double match_distance_in_pixels2;
1374 anbool solved;
1375 double logaccept;
1376
1377 mo->indexid = sp->index->indexid;
1378 mo->healpix = sp->index->healpix;
1379 mo->hpnside = sp->index->hpnside;
1380 mo->wcstan.imagew = sp->field_maxx;
1381 mo->wcstan.imageh = sp->field_maxy;
1382 mo->dimquads = quadfile_dimquads(sp->index->quads);
1383
1384 match_distance_in_pixels2 = square(sp->verify_pix) +
1385 square(sp->index->index_jitter / mo->scale);
1386
1387 logaccept = MIN(sp->logratio_tokeep, sp->logratio_totune);
1388
1389 verify_hit(sp->index->starkd, sp->index->cutnside,
1390 mo, sip, sp->vf, match_distance_in_pixels2,
1391 sp->distractor_ratio, sp->field_maxx, sp->field_maxy,
1392 sp->logratio_bail_threshold, logaccept,
1393 sp->logratio_stoplooking,
1394 sp->distance_from_quad_bonus, fake_match);
1395 mo->nverified = sp->num_verified++;
1396
1397 if (mo->logodds >= sp->best_logodds) {
1398 sp->best_logodds = mo->logodds;
1399 logverb("Got a new best match: logodds %g.\n", mo->logodds);
1400 }
1401
1402 if (mo->logodds >= sp->logratio_totune &&
1403 mo->logodds < sp->logratio_tokeep) {
1404 logverb("Trying to tune up this solution (logodds = %g; %g)...\n",
1405 mo->logodds, exp(mo->logodds));
1406 solver_tweak2(sp, mo, 1, NULL);
1407 logverb("After tuning, logodds = %g (%g)\n",
1408 mo->logodds, exp(mo->logodds));
1409
1410 // Since we tuned up this solution, we can't just accept the
1411 // resulting log-odds at face value.
1412 if (!fake_match) {
1413 verify_hit(sp->index->starkd, sp->index->cutnside,
1414 mo, mo->sip, sp->vf, match_distance_in_pixels2,
1415 sp->distractor_ratio,
1416 sp->field_maxx, sp->field_maxy,
1417 sp->logratio_bail_threshold,
1418 sp->logratio_tokeep,
1419 sp->logratio_stoplooking,
1420 sp->distance_from_quad_bonus,
1421 fake_match);
1422 logverb("Checking tuned result: logodds = %g (%g)\n",
1423 mo->logodds, exp(mo->logodds));
1424 }
1425 }
1426
1427 if (mo->logodds < sp->logratio_toprint)
1428 return FALSE;
1429
1430 update_timeused(sp);
1431 mo->timeused = sp->timeused;
1432
1433 if(log_get_level()>LOG_ERROR) //# Modified by Robert Lancaster for the StellarSolver Internal Library
1434 matchobj_print(mo, log_get_level());
1435
1436 if (mo->logodds < sp->logratio_tokeep)
1437 return FALSE;
1438
1439 logverb("Pixel scale: %g arcsec/pix.\n", mo->scale);
1440 logverb("Parity: %s.\n", (mo->parity ? "neg" : "pos"));
1441
1442 mo->index = sp->index;
1443 mo->index_jitter = sp->index->index_jitter;
1444
1445 if (sp->predistort) {
1446 int i;
1447 double* matchxy;
1448 double* matchxyz;
1449 double* weights;
1450 int N;
1451 int Ngood;
1452 double dx,dy;
1453
1454 // Apply the distortion.
1455 logverb("Applying the distortion pattern and recomputing WCS...\n");
1456
1457 // this includes conflicts and distractors; we won't fill these arrays.
1458 N = mo->nbest;
1459 matchxy = malloc(N * 2 * sizeof(double));
1460 matchxyz = malloc(N * 3 * sizeof(double));
1461 weights = malloc(N * sizeof(double));
1462
1463 Ngood = 0;
1464 for (i=0; i<N; i++) {
1465 double x,y;
1466 if (mo->theta[i] < 0)
1467 continue;
1468 x = starxy_get_x(sp->fieldxy, i);
1469 y = starxy_get_y(sp->fieldxy, i);
1470 sip_pixel_undistortion(sp->predistort, x, y, &dx, &dy);
1471 matchxy[2*Ngood + 0] = dx;
1472 matchxy[2*Ngood + 1] = dy;
1473 memcpy(matchxyz + 3*Ngood, mo->refxyz + 3*mo->theta[i],
1474 3*sizeof(double));
1475 weights[Ngood] = verify_logodds_to_weight(mo->matchodds[i]);
1476
1477 double xx,yy;
1478 Unused anbool ok;
1479 ok = tan_xyzarr2pixelxy(&mo->wcstan, matchxyz+3*Ngood, &xx, &yy);
1480 assert(ok);
1481 logverb("match: ref(%.1f, %.1f) -- img(%.1f, %.1f) --> dist(%.1f, %.1f)\n",
1482 xx, yy, x, y, dx, dy);
1483
1484 Ngood++;
1485 }
1486
1487 if (sp->do_tweak) {
1488 // Compute the SIP solution using the correspondences
1489 // found during verify(), but with the re-distorted positions.
1490 sip_t sip;
1491 memset(&sip, 0, sizeof(sip_t));
1492 sip.a_order = sip.b_order = sp->tweak_aborder;
1493 sip.ap_order = sip.bp_order = sp->tweak_abporder;
1494 sip.wcstan.imagew = solver_field_width(sp);
1495 sip.wcstan.imageh = solver_field_height(sp);
1496 if (sp->set_crpix) {
1497 sip.wcstan.crpix[0] = sp->crpix[0];
1498 sip.wcstan.crpix[1] = sp->crpix[1];
1499 // find matching crval...
1500 sip_pixel_distortion(sp->predistort,
1501 sp->crpix[0], sp->crpix[1], &dx, &dy);
1502 tan_pixelxy2radecarr(&mo->wcstan, dx, dy, sip.wcstan.crval);
1503
1504 } else {
1505 // keep TAN WCS's crval but distort the crpix.
1506 sip.wcstan.crval[0] = mo->wcstan.crval[0];
1507 sip.wcstan.crval[1] = mo->wcstan.crval[1];
1508 sip_pixel_undistortion(sp->predistort,
1509 mo->wcstan.crpix[0], mo->wcstan.crpix[1],
1510 sip.wcstan.crpix+0, sip.wcstan.crpix+1);
1511 }
1512
1513 int doshift = 1;
1514 fit_sip_wcs(matchxyz, matchxy, weights, N, &(sip.wcstan),
1515 sp->tweak_aborder, sp->tweak_abporder, doshift,
1516 &sip);
1517
1518 for (i=0; i<Ngood; i++) {
1519 double xx,yy;
1520 Unused anbool ok;
1521 ok = sip_xyzarr2pixelxy(&sip, matchxyz+3*i, &xx, &yy);
1522 assert(ok);
1523 logverb("match: ref(%.1f, %.1f) -- dist(%.1f, %.1f)\n",
1524 xx, yy, matchxy[2*i+0], matchxy[2*i+1]);
1525 }
1526
1527 } else {
1528 // Compute new TAN WCS...?
1529 fit_tan_wcs_weighted(matchxyz, matchxy, weights, Ngood,
1530 &mo->wcstan, NULL);
1531 if (sp->set_crpix) {
1532 tan_t wcs2;
1533 fit_tan_wcs_move_tangent_point(matchxyz, matchxy, Ngood,
1534 sp->crpix, &mo->wcstan, &wcs2);
1535 fit_tan_wcs_move_tangent_point(matchxyz, matchxy, Ngood,
1536 sp->crpix, &wcs2, &mo->wcstan);
1537 }
1538 }
1539
1540 free(matchxy);
1541 free(matchxyz);
1542 free(weights);
1543
1544 } else if (sp->do_tweak) {
1545 solver_tweak2(sp, mo, sp->tweak_aborder, sip);
1546
1547 } else if (!sip && sp->set_crpix) {
1548 tan_t wcs2;
1549 tan_t wcs3;
1550 fit_tan_wcs_move_tangent_point(mo->quadxyz, mo->quadpix, mo->dimquads,
1551 sp->crpix, &(mo->wcstan), &wcs2);
1552 fit_tan_wcs_move_tangent_point(mo->quadxyz, mo->quadpix, mo->dimquads,
1553 sp->crpix, &wcs2, &wcs3);
1554 memcpy(&(mo->wcstan), &wcs3, sizeof(tan_t));
1555 /*
1556 Good test case:
1557 http://antwrp.gsfc.nasa.gov/apod/image/0912/Geminid2007_pacholka850wp.jpg
1558 solve-field --config backend.cfg Geminid2007_pacholka850wp.xy \
1559 --scale-low 10 --scale-units degwidth -v --no-tweak --continue --new-fits none \
1560 -o 4 --crpix-center --depth 40-45
1561
1562 printf("Original WCS:\n");
1563 tan_print_to(&(mo->wcstan), stdout);
1564 printf("\n");
1565 printf("Moved WCS:\n");
1566 tan_print_to(&wcs2, stdout);
1567 printf("\n");
1568 printf("Moved again WCS:\n");
1569 tan_print_to(&wcs3, stdout);
1570 printf("\n");
1571 */
1572 }
1573
1574 // If the user didn't supply a callback, or if the callback
1575 // returns TRUE, consider it solved.
1576 solved = (!sp->record_match_callback ||
1577 sp->record_match_callback(mo, sp->userdata));
1578
1579 // New best match?
1580 if (!sp->have_best_match || (mo->logodds > sp->best_match.logodds)) {
1581 if (sp->have_best_match)
1582 verify_free_matchobj(&sp->best_match);
1583 memcpy(&sp->best_match, mo, sizeof(MatchObj));
1584 sp->have_best_match = TRUE;
1585 sp->best_index = sp->index;
1586 } else {
1587 verify_free_matchobj(mo);
1588 }
1589
1590 if (solved) {
1591 sp->best_match_solves = TRUE;
1592 return TRUE;
1593 }
1594 return FALSE;
1595 }
1596
solver_new()1597 solver_t* solver_new() {
1598 solver_t* solver = calloc(1, sizeof(solver_t));
1599 solver_set_default_values(solver);
1600 return solver;
1601 }
1602
solver_set_default_values(solver_t * solver)1603 void solver_set_default_values(solver_t* solver) {
1604 memset(solver, 0, sizeof(solver_t));
1605 solver->indexes = pl_new(16);
1606 solver->funits_upper = HUGE_VAL;
1607 solver->logratio_bail_threshold = log(1e-100);
1608 solver->logratio_stoplooking = HUGE_VAL;
1609 solver->logratio_totune = HUGE_VAL;
1610 solver->parity = DEFAULT_PARITY;
1611 solver->codetol = DEFAULT_CODE_TOL;
1612 solver->distractor_ratio = DEFAULT_DISTRACTOR_RATIO;
1613 solver->verify_pix = DEFAULT_VERIFY_PIX;
1614 solver->verify_uniformize = TRUE;
1615 solver->verify_dedup = TRUE;
1616 solver->distance_from_quad_bonus = TRUE;
1617 solver->tweak_aborder = DEFAULT_TWEAK_ABORDER;
1618 solver->tweak_abporder = DEFAULT_TWEAK_ABPORDER;
1619 }
1620
solver_clear_indexes(solver_t * solver)1621 void solver_clear_indexes(solver_t* solver) {
1622 pl_remove_all(solver->indexes);
1623 solver->index = NULL;
1624 }
1625
solver_cleanup(solver_t * solver)1626 void solver_cleanup(solver_t* solver) {
1627 solver_free_field(solver);
1628 pl_free(solver->indexes);
1629 solver->indexes = NULL;
1630 if (solver->have_best_match) {
1631 verify_free_matchobj(&solver->best_match);
1632 solver->have_best_match = FALSE;
1633 }
1634 if (solver->predistort)
1635 sip_free(solver->predistort);
1636 solver->predistort = NULL;
1637 }
1638
solver_free(solver_t * solver)1639 void solver_free(solver_t* solver) {
1640 if (!solver) return;
1641 solver_cleanup(solver);
1642 free(solver);
1643 }
1644