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 #include <assert.h>
7 #include <math.h>
8 #include <string.h>
9 #include <stdint.h>
10
11 #include "verify2.h"
12 #include "permutedsort.h"
13 #include "mathutil.h"
14 #include "keywords.h"
15 #include "log.h"
16 #include "sip-utils.h"
17 #include "healpix.h"
18
19 #define DEBUGVERIFY 1
20 #if DEBUGVERIFY
21 //#define debug(args...) fprintf(stderr, args)
22 #else
23 #define debug(args...)
24 #endif
25
26 #include "fitsioutils.h"
27 #include "errors.h"
28 #include "cairoutils.h"
29
add_gaussian_to_image(double * img,int W,int H,double cx,double cy,double sigma,double scale,double nsigma,int boundary)30 static void add_gaussian_to_image(double* img, int W, int H,
31 double cx, double cy, double sigma,
32 double scale,
33 double nsigma, int boundary) {
34 int x, y;
35 if (boundary == 0) {
36 // truncate.
37 for (y = MAX(0, cy - nsigma*sigma); y <= MIN(H-1, cy + nsigma * sigma); y++) {
38 for (x = MAX(0, cx - nsigma*sigma); x <= MIN(W-1, cx + nsigma * sigma); x++) {
39 img[y*W + x] += scale * exp(-(square(y-cy)+square(x-cx)) / (2.0 * square(sigma)));
40 }
41 }
42
43 } else if (boundary == 1) {
44 // mirror.
45 int mx, my;
46 for (y=MAX(-(H-1), floor(cy - nsigma * sigma));
47 y<=MIN(2*H-1, ceil(cy + nsigma * sigma)); y++) {
48 if (y < 0)
49 my = -1 - y;
50 else if (y >= H)
51 my = 2*H - 1 - y;
52 else
53 my = y;
54 for (x=MAX(-(W-1), floor(cx - nsigma * sigma));
55 x<=MIN(2*W-1, ceil(cx + nsigma * sigma)); x++) {
56 if (x < 0)
57 mx = -1 - x;
58 else if (x >= W)
59 mx = 2*W - 1 - x;
60 else
61 mx = x;
62 img[my*W + mx] += scale * exp(-(square(y-cy)+square(x-cx)) / (2.0 * square(sigma)));
63 }
64 }
65 }
66 }
67
verify_get_all_matches(const double * refxys,int NR,const double * testxys,const double * testsigma2s,int NT,double effective_area,double distractors,double nsigma,double limit,il *** p_reflist,dl *** p_problist)68 void verify_get_all_matches(const double* refxys, int NR,
69 const double* testxys, const double* testsigma2s, int NT,
70 double effective_area,
71 double distractors,
72 double nsigma,
73 double limit,
74 il*** p_reflist,
75 dl*** p_problist) {
76 double* refcopy;
77 kdtree_t* rtree;
78 int Nleaf = 10;
79 int i,j;
80 double logd;
81 double logbg;
82 double loglimit;
83
84 il** reflist;
85 dl** problist;
86
87 reflist = calloc(NT, sizeof(il*));
88 problist = calloc(NT, sizeof(dl*));
89
90 // Build a tree out of the index stars in pixel space...
91 // kdtree scrambles the data array so make a copy first.
92 refcopy = malloc(2 * NR * sizeof(double));
93 memcpy(refcopy, refxys, 2 * NR * sizeof(double));
94 rtree = kdtree_build(NULL, refcopy, NR, 2, Nleaf, KDTT_DOUBLE, KD_BUILD_SPLIT);
95
96 logbg = log(1.0 / effective_area);
97 logd = log(distractors / effective_area);
98 loglimit = log(distractors / effective_area * limit);
99
100 for (i=0; i<NT; i++) {
101 const double* testxy;
102 double sig2;
103 kdtree_qres_t* res;
104
105 testxy = testxys + 2*i;
106 sig2 = testsigma2s[i];
107
108 logverb("\n");
109 logverb("test star %i: (%.1f,%.1f), sigma: %.1f\n", i, testxy[0], testxy[1], sqrt(sig2));
110
111 // find all ref stars within nsigma.
112 res = kdtree_rangesearch_options(rtree, testxy, sig2*nsigma*nsigma,
113 KD_OPTIONS_SORT_DISTS | KD_OPTIONS_SMALL_RADIUS);
114
115 if (res->nres == 0) {
116 kdtree_free_query(res);
117 continue;
118 }
119
120 reflist[i] = il_new(4);
121 problist[i] = dl_new(4);
122
123 for (j=0; j<res->nres; j++) {
124 double d2;
125 int refi;
126 double loggmax, logfg;
127
128 d2 = res->sdists[j];
129 refi = res->inds[j];
130
131 // peak value of the Gaussian
132 loggmax = log((1.0 - distractors) / (2.0 * M_PI * sig2 * NR));
133 // value of the Gaussian
134 logfg = loggmax - d2 / (2.0 * sig2);
135
136 if (logfg < loglimit)
137 continue;
138
139 logverb(" ref star %i, dist %.2f, sigmas: %.3f, logfg: %.1f (%.1f above distractor, %.1f above bg, %.1f above keep-limit)\n",
140 refi, sqrt(d2), sqrt(d2 / sig2), logfg, logfg - logd, logfg - logbg, logfg - loglimit);
141
142 il_append(reflist[i], refi);
143 dl_append(problist[i], logfg);
144 }
145
146 kdtree_free_query(res);
147 }
148 kdtree_free(rtree);
149 free(refcopy);
150
151 *p_reflist = reflist;
152 *p_problist = problist;
153 }
154
155
156 /*
157 unsigned char* img = malloc(4*W*H);
158
159 // draw images of index and field densities.
160 double* idensity = calloc(W * H, sizeof(double));
161 double* fdensity = calloc(W * H, sizeof(double));
162
163 double iscale = 2. * sqrt((double)(W * H) / (NI * M_PI));
164 double fscale = 2. * sqrt((double)(W * H) / (NF * M_PI));
165 logverb("NI = %i; iscale = %g\n", NI, iscale);
166 logverb("NF = %i; fscale = %g\n", NF, fscale);
167 logverb("computing density images...\n");
168 for (i=0; i<NI; i++)
169 add_gaussian_to_image(idensity, W, H, indexpix[i*2 + 0], indexpix[i*2 + 1], iscale, 1.0, 3.0, 1);
170 for (i=0; i<NF; i++)
171 add_gaussian_to_image(fdensity, W, H, starxy_getx(vf->field, i), starxy_gety(vf->field, i), fscale, 1.0, 3.0, 1);
172
173 double idmax=0, fdmax=0;
174 for (i=0; i<(W*H); i++) {
175 idmax = MAX(idmax, idensity[i]);
176 fdmax = MAX(fdmax, fdensity[i]);
177 }
178 for (i=0; i<(W*H); i++) {
179 unsigned char val = 255.5 * idensity[i] / idmax;
180 img[i*4+0] = val;
181 img[i*4+1] = val;
182 img[i*4+2] = val;
183 img[i*4+3] = 255;
184 }
185 cairoutils_write_png("idensity.png", img, W, H);
186 for (i=0; i<(W*H); i++) {
187 unsigned char val = 255.5 * fdensity[i] / fdmax;
188 img[i*4+0] = val;
189 img[i*4+1] = val;
190 img[i*4+2] = val;
191 img[i*4+3] = 255;
192 }
193 cairoutils_write_png("fdensity.png", img, W, H);
194
195 free(idensity);
196 free(fdensity);
197 */
198
199
200 /*
201 // index star weights.
202 double* iweights = malloc(NI * sizeof(double));
203 for (i=0; i<NI; i++)
204 iweights[i] = 1.0;
205 // don't count index stars that are part of the matched quad.
206 for (i=0; i<NI; i++)
207 for (j=0; j<dimquads; j++)
208 if (starids[i] == mo->star[j])
209 iweights[i] = 0.0;
210
211 double ranksigma = 20.0;
212 int Nrankprobs = (int)(ranksigma * 5);
213 double* rankprobs = malloc(Nrankprobs * sizeof(double));
214 for (i=0; i<Nrankprobs; i++)
215 rankprobs[i] = exp(-(double)(i*i) / (2. * ranksigma * ranksigma));
216
217 double qc[2], qr2;
218
219 get_quad_center(vf, mo, qc, &qr2);
220 */
221
222 /*
223 // create the probability distribution map for this field star.
224 double* pmap = malloc(W * H * sizeof(double));
225 // background rate...
226 for (j=0; j<(W*H); j++)
227 pmap[j] = distractors / (fieldW*fieldH);
228
229 double normrankprob;
230 int dr;
231 normrankprob = 0.0;
232 for (j=0; j<NI; j++) {
233 dr = abs(i - j);
234 if (dr >= Nrankprobs)
235 continue;
236 normrankprob += rankprobs[dr];
237 }
238 for (j=0; j<NI; j++) {
239 double r2, sig2, sig;
240
241 dr = abs(i - j);
242 if (dr >= Nrankprobs)
243 continue;
244
245 r2 = distsq(indexpix + j*2, qc, 2);
246 sig2 = get_sigma2_at_radius(verify_pix2, r2, qr2);
247 sig = sqrt(sig2);
248
249 for (y = MAX(0, indexpix[j*2 + 1] - 5.0*sig);
250 y <= MIN(H-1, indexpix[j*2 + 1] + 5.0 * sig); y++) {
251 for (x = MAX(0, indexpix[j*2 + 0] - 5.0*sig);
252 x <= MIN(W-1, indexpix[j*2 + 0] + 5.0 * sig); x++) {
253 pmap[y*W + x] += iweights[j] * (rankprobs[dr] / normrankprob) *
254 exp(-(square(indexpix[j*2+0]-x) + square(indexpix[j*2+1]-y)) / (2.0 * sig2)) *
255 (1.0 / (2.0 * M_PI * sig2)) * (1.0 - distractors);
256 }
257 }
258 }
259
260 double maxval = 0.0;
261 for (j=0; j<(W*H); j++)
262 maxval = MAX(maxval, pmap[j]);
263
264 double psum = 0.0;
265 for (j=0; j<(W*H); j++)
266 psum += pmap[j];
267 logverb("Probability sum: %g\n", psum);
268
269 char* fn;
270
271 printf("Writing probability image %i\n", i);
272 printf("maxval = %g\n", maxval);
273
274 for (j=0; j<(W*H); j++) {
275 unsigned char val = (unsigned char)MAX(0, MIN(255, 255.5 * sqrt(pmap[j] / maxval)));
276 img[4*j + 0] = val;
277 img[4*j + 1] = val;
278 img[4*j + 2] = val;
279 img[4*j + 3] = 255;
280 }
281
282 free(pmap);
283
284 cairo_t* cairo;
285 cairo_surface_t* target;
286
287 target = cairo_image_surface_create_for_data(img, CAIRO_FORMAT_ARGB32, W, H, W*4);
288 cairo = cairo_create(target);
289 cairo_set_line_join(cairo, CAIRO_LINE_JOIN_ROUND);
290 cairo_set_antialias(cairo, CAIRO_ANTIALIAS_GRAY);
291
292 // draw the quad.
293 cairo_set_source_rgba(cairo, 0.0, 0.0, 1.0, 0.5);
294
295 double starxy[DQMAX*2];
296 double angles[DQMAX];
297 int perm[DQMAX];
298 int DQ = dimquads;
299
300 for (k=0; k<DQ; k++)
301 starxy_get(vf->field, mo->field[k], starxy + k*2);
302 for (k=0; k<DQ; k++)
303 angles[k] = atan2(starxy[k*2 + 1] - qc[1], starxy[k*2 + 0] - qc[0]);
304 permutation_init(perm, DQ);
305 permuted_sort(angles, sizeof(double), compare_doubles_asc, perm, DQ);
306 for (k=0; k<DQ; k++) {
307 double px, py;
308 px = starxy[perm[k]*2+0];
309 py = starxy[perm[k]*2+1];
310 if (k == 0)
311 cairo_move_to(cairo, px, py);
312 else
313 cairo_line_to(cairo, px, py);
314 }
315 cairo_close_path(cairo);
316 cairo_stroke(cairo);
317
318 // draw the source point.
319 cairo_set_source_rgb(cairo, 1.0, 0.0, 0.0);
320 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_CROSSHAIR, fxy[0], fxy[1], 3.0);
321 cairo_stroke(cairo);
322
323 */
324
325
326 /*
327 // find index stars within n sigma...
328 double nsigma = 5.0;
329 double r2 = sigma2s[i] * nsigma*nsigma;
330 int opts = KD_OPTIONS_COMPUTE_DISTS;
331
332 kdtree_qres_t* res = kdtree_rangesearch_options(itree, fxy, r2, opts);
333
334 logverb("found %i index stars within range.\n", res->nres);
335
336 if (res->nres == 0)
337 goto nextstar;
338
339 double bestprob = 0.0;
340 int bestarg = -1;
341
342 for (j=0; j<res->nres; j++) {
343 int ii = res->inds[j];
344 dr = abs(i - ii);
345 logverb(" index star %i: rank diff %i\n", ii, i-ii);
346 if (dr >= Nrankprobs)
347 continue;
348 double prob = iweights[ii] * rankprobs[dr] *
349 exp(-res->sdists[j] / (2.0 * sigma2s[i]));
350 logverb(" -> prob %g\n", prob);
351 if (prob > bestprob) {
352 bestprob = prob;
353 bestarg = j;
354 }
355 }
356
357 if (bestarg == -1) {
358 // FIXME -- distractor?
359 goto nextstar;
360 }
361
362 int besti = res->inds[bestarg];
363 double lostweight = bestprob / iweights[besti];
364 // exp(-res->sdists[bestarg] / (2.0 * sigma2s[i]));
365 iweights[besti] = MAX(0.0, iweights[besti] - lostweight);
366 logverb("matched index star %i: dropping weight by %g; new weight %g.\n", besti, lostweight, iweights[besti]);
367
368 bestprob *= (1.0 / (2.0 * M_PI * sigma2s[i])) / normrankprob * (1.0 - distractors);
369 logverb("bestprob %g, vs background %g\n", bestprob, (1.0 / (double)(W * H)));
370
371 kdtree_free_query(res);
372
373 cairo_set_source_rgb(cairo, 0.0, 1.0, 0.0);
374 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_CIRCLE, indexpix[2*besti+0], indexpix[2*besti+1], 3.0);
375 cairo_stroke(cairo);
376
377
378 nextstar:
379 cairoutils_argb32_to_rgba(img, W, H);
380 cairo_surface_destroy(target);
381 cairo_destroy(cairo);
382 asprintf(&fn, "logprob-%04i.png", i);
383 cairoutils_write_png(fn, img, W, H);
384 free(img);
385
386
387 if (i >= 9)
388 break;
389 }
390
391 kdtree_free(itree);
392 free(ipixcopy);
393
394 */
395
396
397
398
399
400
401
402