1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5 #include <stdio.h>
6 #include <assert.h>
7
8 #include "os-features.h"
9 #include "bl.h"
10 #include "fit-wcs.h"
11 #include "sip.h"
12 #include "sip_qfits.h"
13 #include "sip-utils.h"
14 #include "scamp.h"
15 #include "log.h"
16 #include "errors.h"
17 #include "tweak.h"
18 #include "matchfile.h"
19 #include "matchobj.h"
20 #include "boilerplate.h"
21 #include "xylist.h"
22 #include "rdlist.h"
23 #include "mathutil.h"
24 #include "verify.h"
25 #include "fitsioutils.h"
26
27
28 // Tweak debug plots?
29 #define TWEAK_DEBUG_PLOTS 0
30 #if TWEAK_DEBUG_PLOTS
31 #include "plotstuff.h"
32 #include "plotimage.h"
33 #include "cairoutils.h"
34
makeplot(const char * plotfn,char * bgimgfn,int W,int H,int Nfield,double * fieldpix,double * fieldsigma2s,int Nindex,double * indexpix,int besti,int * theta,double * crpix,int * testperm,double * qc)35 static void makeplot(const char* plotfn, char* bgimgfn, int W, int H,
36 int Nfield, double* fieldpix, double* fieldsigma2s,
37 int Nindex, double* indexpix, int besti, int* theta,
38 double* crpix, int* testperm,
39 double * qc) {
40 int i;
41 plot_args_t pargs;
42 plotimage_t* img;
43 cairo_t* cairo;
44 int ti;
45 logmsg("Creating plot %s\n", plotfn);
46 plotstuff_init(&pargs);
47 pargs.outformat = PLOTSTUFF_FORMAT_PNG;
48 pargs.outfn = plotfn;
49 pargs.fontsize = 12;
50 if (bgimgfn) {
51 img = plotstuff_get_config(&pargs, "image");
52 img->format = PLOTSTUFF_FORMAT_JPG;
53 plot_image_set_filename(img, bgimgfn);
54 plot_image_setsize(&pargs, img);
55 plotstuff_run_command(&pargs, "image");
56 } else {
57 float rgba[4] = {0, 0, 0.1, 1.0};
58 plotstuff_set_size(&pargs, W, H);
59 //plotstuff_init2(&pargs);
60 plotstuff_set_rgba(&pargs, rgba);
61 plotstuff_run_command(&pargs, "fill");
62 }
63 cairo = pargs.cairo;
64 // red circles around every field star.
65 cairo_set_color(cairo, "gray");
66 for (i=0; i<Nfield; i++) {
67 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_CIRCLE,
68 fieldpix[2*i+0], fieldpix[2*i+1],
69 2.0 * sqrt(fieldsigma2s[i]));
70 cairo_stroke(cairo);
71 }
72 // green crosshairs at every index star.
73 cairo_set_color(cairo, "green");
74 for (i=0; i<Nindex; i++) {
75 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_XCROSSHAIR,
76 indexpix[2*i+0], indexpix[2*i+1], 3);
77 cairo_stroke(cairo);
78 }
79 // thick white circles for corresponding field stars.
80 cairo_set_line_width(cairo, 2);
81 for (ti=0; ti<=besti; ti++) {
82 if (testperm)
83 i = testperm[ti];
84 else
85 i = ti;
86 //printf("field %i -> index %i\n", i, theta[i]);
87 if (theta[i] < 0)
88 continue;
89 cairo_set_color(cairo, "white");
90 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_CIRCLE,
91 fieldpix[2*i+0], fieldpix[2*i+1],
92 2.0 * sqrt(fieldsigma2s[i]));
93 cairo_stroke(cairo);
94 // thick cyan crosshairs for corresponding index stars.
95 cairo_set_color(cairo, "cyan");
96 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_XCROSSHAIR,
97 indexpix[2*theta[i]+0],
98 indexpix[2*theta[i]+1],
99 3);
100 cairo_stroke(cairo);
101 }
102
103 cairo_set_line_width(cairo, 2);
104
105 //for (i=0; i<=besti; i++) {
106 for (ti=0; ti<Nfield; ti++) {
107 anbool mark = TRUE;
108 if (testperm)
109 i = testperm[ti];
110 else
111 i = ti;
112 switch (theta[i]) {
113 case THETA_DISTRACTOR:
114 cairo_set_color(cairo, "red");
115 break;
116 case THETA_CONFLICT:
117 cairo_set_color(cairo, "yellow");
118 break;
119 case THETA_FILTERED:
120 cairo_set_color(cairo, "orange");
121 break;
122 default:
123 if (theta[i] < 0) {
124 cairo_set_color(cairo, "gray");
125 } else {
126 cairo_set_color(cairo, "white");
127 }
128 mark = FALSE;
129 }
130
131 if (mark) {
132 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_CIRCLE,
133 fieldpix[2*i+0], fieldpix[2*i+1],
134 2.0 * sqrt(fieldsigma2s[i]));
135 cairo_stroke(cairo);
136 }
137
138 if (ti <= MAX(besti, 10)) {
139 char label[32];
140 sprintf(label, "%i", i);
141 plotstuff_text_xy(&pargs, fieldpix[2*i+0], fieldpix[2*i+1], label);
142 }
143 if (i == besti) {
144 cairo_set_line_width(cairo, 1);
145 }
146 }
147
148
149 if (crpix) {
150 cairo_set_color(cairo, "yellow");
151 cairo_set_line_width(cairo, 4);
152 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_CROSSHAIR,
153 crpix[0], crpix[1], 10);
154 cairo_stroke(cairo);
155 }
156
157 if (qc) {
158 cairo_set_color(cairo, "skyblue");
159 cairo_set_line_width(cairo, 4);
160 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_CROSSHAIR,
161 qc[0], qc[1], 10);
162 cairo_stroke(cairo);
163 }
164
165 plotstuff_output(&pargs);
166 logmsg("Wrote plot %s\n", plotfn);
167 }
168
tdebugfn(const char * name)169 static char* tdebugfn(const char* name) {
170 static char fn[256];
171 static int plotnum = 0;
172 sprintf(fn, "tweak-%03i-%s.png", plotnum, name);
173 plotnum++;
174 return fn;
175 }
176
177 #define TWEAK_DEBUG_PLOT(name, W, H, Nfield, fieldxy, fieldsig2, \
178 Nindex, indexxy, besti, theta, crpix, testperm, qc) \
179 makeplot(tdebugfn(name), NULL, W, H, Nfield, fieldxy, fieldsig2, \
180 Nindex, indexxy, besti, theta, crpix, testperm, qc);
181
182
183 #else
184
185 #define TWEAK_DEBUG_PLOT(name, W, H, Nfield, fieldxy, fieldsig2, \
186 Nindex, indexxy, besti, theta, crpix, testperm, qc) \
187 do{}while(0)
188
189 #endif
190
191
192
193
194
tweak2(const double * fieldxy,int Nfield,double fieldjitter,int W,int H,const double * indexradec,int Nindex,double indexjitter,const double * quadcenter,double quadR2,double distractors,double logodds_bail,int sip_order,int sip_invorder,const sip_t * startwcs,sip_t * destwcs,int ** newtheta,double ** newodds,double * crpix,double * p_logodds,int * p_besti,int * testperm,int startorder)195 sip_t* tweak2(const double* fieldxy, int Nfield,
196 double fieldjitter,
197 int W, int H,
198 const double* indexradec, int Nindex,
199 double indexjitter,
200 const double* quadcenter, double quadR2,
201 double distractors,
202 double logodds_bail,
203 int sip_order,
204 int sip_invorder,
205 const sip_t* startwcs,
206 sip_t* destwcs,
207 int** newtheta, double** newodds,
208 double* crpix,
209 double* p_logodds,
210 int* p_besti,
211 int* testperm,
212 int startorder) {
213 int order;
214 sip_t* sipout;
215 int* indexin;
216 double* indexpix;
217 double* fieldsigma2s;
218 double* weights;
219 double* matchxyz;
220 double* matchxy;
221 int i, Nin=0;
222 double logodds = 0;
223 int besti = -1;
224 int* theta = NULL;
225 double* odds = NULL;
226 int* refperm = NULL;
227 double qc[2];
228
229 memcpy(qc, quadcenter, 2*sizeof(double));
230
231 if (destwcs)
232 sipout = destwcs;
233 else
234 sipout = sip_create();
235
236 indexin = malloc(Nindex * sizeof(int));
237 indexpix = malloc(2 * Nindex * sizeof(double));
238 fieldsigma2s = malloc(Nfield * sizeof(double));
239 weights = malloc(Nfield * sizeof(double));
240 matchxyz = malloc(Nfield * 3 * sizeof(double));
241 matchxy = malloc(Nfield * 2 * sizeof(double));
242
243 // FIXME --- hmmm, how do the annealing steps and iterating up to
244 // higher orders interact?
245
246 assert(startwcs);
247 memcpy(sipout, startwcs, sizeof(sip_t));
248
249 logverb("tweak2: starting orders %i, %i\n", sipout->a_order, sipout->ap_order);
250
251 if (!sipout->wcstan.imagew)
252 sipout->wcstan.imagew = W;
253 if (!sipout->wcstan.imageh)
254 sipout->wcstan.imageh = H;
255
256 logverb("Tweak2: starting from WCS:\n");
257 if (log_get_level() >= LOG_VERB)
258 sip_print_to(sipout, stdout);
259
260 for (order=startorder; order <= sip_order; order++) {
261 int step;
262 int STEPS = 100;
263 // variance growth rate wrt radius.
264 double gamma = 1.0;
265 //logverb("Starting tweak2 order=%i\n", order);
266
267 for (step=0; step<STEPS; step++) {
268 double iscale;
269 double ijitter;
270 double ra, dec;
271 double R2;
272 int Nmatch;
273 int nmatch, nconf, ndist;
274 double pix2;
275 double totalweight;
276
277 // clean up from last round (we do it here so that they're
278 // valid when we leave the loop)
279 free(theta);
280 free(odds);
281 free(refperm);
282
283 // Anneal
284 gamma = pow(0.9, step);
285 if (step == STEPS-1)
286 gamma = 0.0;
287 logverb("Annealing: order %i, step %i, gamma = %g\n", order, step, gamma);
288
289 debug("Using input WCS:\n");
290 if (log_get_level() > LOG_VERB)
291 sip_print_to(sipout, stdout);
292
293 // Project reference sources into pixel space; keep the ones inside image bounds.
294 Nin = 0;
295 for (i=0; i<Nindex; i++) {
296 anbool ok;
297 double x,y;
298 ra = indexradec[2*i + 0];
299 dec = indexradec[2*i + 1];
300 ok = sip_radec2pixelxy(sipout, ra, dec, &x, &y);
301 if (!ok)
302 continue;
303 if (!sip_pixel_is_inside_image(sipout, x, y))
304 continue;
305 indexpix[Nin*2+0] = x;
306 indexpix[Nin*2+1] = y;
307 indexin[Nin] = i;
308 Nin++;
309 }
310 logverb("%i reference sources within the image.\n", Nin);
311 //logverb("CRPIX is (%g,%g)\n", sip.wcstan.crpix[0], sip.wcstan.crpix[1]);
312
313 if (Nin == 0) {
314 sip_free(sipout);
315 free(matchxy);
316 free(matchxyz);
317 free(weights);
318 free(fieldsigma2s);
319 free(indexpix);
320 free(indexin);
321 return NULL;
322 }
323
324 iscale = sip_pixel_scale(sipout);
325 ijitter = indexjitter / iscale;
326 //logverb("With pixel scale of %g arcsec/pixel, index adds jitter of %g pix.\n", iscale, ijitter);
327
328 /* CHECK
329 for (i=0; i<Nin; i++) {
330 double x,y;
331 int ii = indexin[i];
332 sip_radec2pixelxy(sipout, indexradec[2*ii+0], indexradec[2*ii+1], &x, &y);
333 logverb("indexin[%i]=%i; (%.1f,%.1f) -- (%.1f,%.1f)\n",
334 i, ii, indexpix[i*2+0], indexpix[i*2+1], x, y);
335 }
336 */
337
338 for (i=0; i<Nfield; i++) {
339 R2 = distsq(qc, fieldxy + 2*i, 2);
340 fieldsigma2s[i] = (square(fieldjitter) + square(ijitter)) * (1.0 + gamma * R2/quadR2);
341 }
342
343 if (order == 1 && step == 0 && TWEAK_DEBUG_PLOTS) {
344 TWEAK_DEBUG_PLOT("init", W, H, Nfield, fieldxy, fieldsigma2s,
345 Nin, indexpix, *p_besti, *newtheta,
346 sipout->wcstan.crpix, testperm, qc);
347 }
348
349 /*
350 logodds = verify_star_lists(indexpix, Nin,
351 fieldxy, fieldsigma2s, Nfield,
352 W*H, distractors,
353 logodds_bail, HUGE_VAL,
354 &besti, &odds, &theta, NULL,
355 &testperm);
356 */
357
358 pix2 = square(fieldjitter);
359 logodds = verify_star_lists_ror(indexpix, Nin,
360 fieldxy, fieldsigma2s, Nfield,
361 pix2, gamma, qc, quadR2,
362 W, H, distractors,
363 logodds_bail, HUGE_VAL,
364 &besti, &odds, &theta, NULL,
365 &testperm, &refperm);
366
367 logverb("Logodds: %g\n", logodds);
368 verify_count_hits(theta, besti, &nmatch, &nconf, &ndist);
369 logverb("%i matches, %i distractors, %i conflicts (at best log-odds); %i field sources, %i index sources\n", nmatch, ndist, nconf, Nfield, Nin);
370 verify_count_hits(theta, Nfield-1, &nmatch, &nconf, &ndist);
371 logverb("%i matches, %i distractors, %i conflicts (all sources)\n", nmatch, ndist, nconf);
372 if (log_get_level() >= LOG_VERB) {
373 matchobj_log_hit_miss(theta, testperm, besti+1, Nfield, LOG_VERB, "Hit/miss: ");
374 }
375
376 /*
377 logverb("\nAfter verify():\n");
378 for (i=0; i<Nin; i++) {
379 double x,y;
380 int ii = indexin[refperm[i]];
381 sip_radec2pixelxy(sipout, indexradec[2*ii+0], indexradec[2*ii+1], &x, &y);
382 logverb("indexin[%i]=%i; (%.1f,%.1f) -- (%.1f,%.1f)\n",
383 i, ii, indexpix[i*2+0], indexpix[i*2+1], x, y);
384 }
385 */
386
387 if (TWEAK_DEBUG_PLOTS) {
388 char name[32];
389 sprintf(name, "o%is%02ipre", order, step);
390 TWEAK_DEBUG_PLOT(name, W, H, Nfield, fieldxy, fieldsigma2s,
391 Nin, indexpix, besti, theta,
392 sipout->wcstan.crpix, testperm, qc);
393 }
394
395 Nmatch = 0;
396 debug("Weights:");
397 for (i=0; i<Nfield; i++) {
398 double ra,dec;
399 if (theta[i] < 0)
400 continue;
401 assert(theta[i] < Nin);
402 int ii = indexin[refperm[theta[i]]];
403 assert(ii < Nindex);
404 assert(ii >= 0);
405
406 ra = indexradec[ii*2+0];
407 dec = indexradec[ii*2+1];
408 radecdeg2xyzarr(ra, dec, matchxyz + Nmatch*3);
409 memcpy(matchxy + Nmatch*2, fieldxy + i*2, 2*sizeof(double));
410 weights[Nmatch] = verify_logodds_to_weight(odds[i]);
411 debug(" %.2f", weights[Nmatch]);
412 Nmatch++;
413
414 /*
415 logverb("match img (%.1f,%.1f) -- ref (%.1f, %.1f), odds %g, wt %.3f\n",
416 fieldxy[i*2+0], fieldxy[i*2+1],
417 indexpix[theta[i]*2+0], indexpix[theta[i]*2+1],
418 odds[i],
419 weights[Nmatch-1]);
420 double xx,yy;
421 sip_radec2pixelxy(sipout, ra, dec, &xx, &yy);
422 logverb("check: (%.1f, %.1f)\n", xx, yy);
423 */
424 }
425 debug("\n");
426
427 if (Nmatch < 2) {
428 logverb("No matches -- aborting tweak attempt\n");
429 free(theta);
430 sip_free(sipout);
431 free(matchxy);
432 free(matchxyz);
433 free(weights);
434 free(fieldsigma2s);
435 free(indexpix);
436 free(indexin);
437 return NULL;
438 }
439
440 // Update the "quad center" to be the weighted average matched star posn.
441 qc[0] = qc[1] = 0.0;
442 totalweight = 0.0;
443 for (i=0; i<Nmatch; i++) {
444 qc[0] += (weights[i] * matchxy[2*i+0]);
445 qc[1] += (weights[i] * matchxy[2*i+1]);
446 totalweight += weights[i];
447 }
448 qc[0] /= totalweight;
449 qc[1] /= totalweight;
450 logverb("Moved quad center to (%.1f, %.1f)\n", qc[0], qc[1]);
451
452 //
453 sipout->a_order = sipout->b_order = order;
454 sipout->ap_order = sipout->bp_order = sip_invorder;
455 logverb("tweak2: setting orders %i, %i\n", sipout->a_order, sipout->ap_order);
456
457 if (crpix) {
458 tan_t temptan;
459 logverb("Moving tangent point to given CRPIX (%g,%g)\n", crpix[0], crpix[1]);
460 fit_tan_wcs_move_tangent_point_weighted(matchxyz, matchxy, weights, Nmatch,
461 crpix, &sipout->wcstan, &temptan);
462 fit_tan_wcs_move_tangent_point_weighted(matchxyz, matchxy, weights, Nmatch,
463 crpix, &temptan, &sipout->wcstan);
464 }
465
466 int doshift = 1;
467 fit_sip_wcs(matchxyz, matchxy, weights, Nmatch,
468 &(sipout->wcstan), order, sip_invorder,
469 doshift, sipout);
470
471 debug("Got SIP:\n");
472 if (log_get_level() > LOG_VERB)
473 sip_print_to(sipout, stdout);
474 sipout->wcstan.imagew = W;
475 sipout->wcstan.imageh = H;
476 }
477 }
478
479 //logverb("Final logodds: %g\n", logodds);
480
481 // Now, recompute final logodds after turning 'gamma' on again (?)
482 // FIXME -- this counts the quad stars in the logodds...
483 {
484 double gamma = 1.0;
485 double iscale;
486 double ijitter;
487 double ra, dec;
488 double R2;
489 int nmatch, nconf, ndist;
490 double pix2;
491
492 free(theta);
493 free(odds);
494 free(refperm);
495 gamma = 1.0;
496 // Project reference sources into pixel space; keep the ones inside image bounds.
497 Nin = 0;
498 for (i=0; i<Nindex; i++) {
499 anbool ok;
500 double x,y;
501 ra = indexradec[2*i + 0];
502 dec = indexradec[2*i + 1];
503 ok = sip_radec2pixelxy(sipout, ra, dec, &x, &y);
504 if (!ok)
505 continue;
506 if (!sip_pixel_is_inside_image(sipout, x, y))
507 continue;
508 indexpix[Nin*2+0] = x;
509 indexpix[Nin*2+1] = y;
510 indexin[Nin] = i;
511 Nin++;
512 }
513 logverb("%i reference sources within the image.\n", Nin);
514
515 iscale = sip_pixel_scale(sipout);
516 ijitter = indexjitter / iscale;
517 for (i=0; i<Nfield; i++) {
518 R2 = distsq(qc, fieldxy + 2*i, 2);
519 fieldsigma2s[i] = (square(fieldjitter) + square(ijitter)) * (1.0 + gamma * R2/quadR2);
520 }
521
522 pix2 = square(fieldjitter);
523 logodds = verify_star_lists_ror(indexpix, Nin,
524 fieldxy, fieldsigma2s, Nfield,
525 pix2, gamma, qc, quadR2,
526 W, H, distractors,
527 logodds_bail, HUGE_VAL,
528 &besti, &odds, &theta, NULL,
529 &testperm, &refperm);
530 logverb("Logodds: %g\n", logodds);
531 verify_count_hits(theta, besti, &nmatch, &nconf, &ndist);
532 logverb("%i matches, %i distractors, %i conflicts (at best log-odds); %i field sources, %i index sources\n", nmatch, ndist, nconf, Nfield, Nin);
533 verify_count_hits(theta, Nfield-1, &nmatch, &nconf, &ndist);
534 logverb("%i matches, %i distractors, %i conflicts (all sources)\n", nmatch, ndist, nconf);
535 if (log_get_level() >= LOG_VERB) {
536 matchobj_log_hit_miss(theta, testperm, besti+1, Nfield, LOG_VERB,
537 "Hit/miss: ");
538 }
539
540 if (TWEAK_DEBUG_PLOTS) {
541 TWEAK_DEBUG_PLOT("final", W, H, Nfield, fieldxy, fieldsigma2s,
542 Nin, indexpix, besti, theta,
543 sipout->wcstan.crpix, testperm, qc);
544 }
545 }
546
547
548 if (newtheta) {
549 // undo the "indexpix" inside-image-bounds cut.
550 (*newtheta) = malloc(Nfield * sizeof(int));
551 for (i=0; i<Nfield; i++) {
552 int nt;
553 if (theta[i] < 0)
554 nt = theta[i];
555 else
556 nt = indexin[refperm[theta[i]]];
557 (*newtheta)[i] = nt;
558 }
559 }
560 free(theta);
561 free(refperm);
562
563 if (newodds)
564 *newodds = odds;
565 else
566 free(odds);
567
568 logverb("Tweak2: final WCS:\n");
569 if (log_get_level() >= LOG_VERB)
570 sip_print_to(sipout, stdout);
571
572 if (p_logodds)
573 *p_logodds = logodds;
574 if (p_besti)
575 *p_besti = besti;
576
577 free(indexin);
578 free(indexpix);
579 free(fieldsigma2s);
580 free(weights);
581 free(matchxyz);
582 free(matchxy);
583
584 return sipout;
585 }
586
587