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 <stdio.h>
7 #include <stdlib.h>
8 #include <assert.h>
9
10 #include "starkd.h"
11 #include "kdtree.h"
12 #include "kdtree_fits_io.h"
13 #include "starutil.h"
14 #include "fitsbin.h"
15 #include "fitstable.h"
16 #include "errors.h"
17 #include "tic.h"
18 #include "log.h"
19 #include "ioutils.h"
20 #include "fitsioutils.h"
21
22 #ifdef _MSC_VER //# Modified by Robert Lancaster for the StellarSolver Internal Library
23 struct timeval {
24 long tv_sec; /* seconds */
25 long tv_usec; /* and microseconds */
26 };
27 #endif
28
startree_alloc()29 static startree_t* startree_alloc() {
30 startree_t* s = calloc(1, sizeof(startree_t));
31 if (!s) {
32 fprintf(stderr, "Failed to allocate a star kdtree struct.\n");
33 return NULL;
34 }
35 return s;
36 }
37
startree_get_tagalong_column_names(startree_t * s,sl * lst)38 sl* startree_get_tagalong_column_names(startree_t* s, sl* lst) {
39 if (!startree_has_tagalong(s))
40 return NULL;
41 return fitstable_get_fits_column_names(startree_get_tagalong(s), lst);
42 }
43
startree_get_tagalong_N_columns(startree_t * s)44 int startree_get_tagalong_N_columns(startree_t* s) {
45 if (!startree_has_tagalong(s))
46 return 0;
47 return fitstable_get_N_fits_columns(startree_get_tagalong(s));
48 }
49
50 /**
51 Returns the name of the 'i'th column in the tagalong table.
52 */
startree_get_tagalong_column_name(startree_t * s,int i)53 const char* startree_get_tagalong_column_name(startree_t* s, int i) {
54 if (!startree_has_tagalong(s))
55 return NULL;
56 return fitstable_get_fits_column_name(startree_get_tagalong(s), i);
57 }
58
startree_get_tagalong_column_fits_type(startree_t * s,int i)59 tfits_type startree_get_tagalong_column_fits_type(startree_t* s, int i) {
60 if (!startree_has_tagalong(s))
61 return TFITS_BIN_TYPE_UNKNOWN;
62 return fitstable_get_fits_column_type(startree_get_tagalong(s), i);
63 }
64
startree_get_tagalong_column_array_size(startree_t * s,int i)65 int startree_get_tagalong_column_array_size(startree_t* s, int i) {
66 if (!startree_has_tagalong(s))
67 return -1;
68 return fitstable_get_fits_column_array_size(startree_get_tagalong(s), i);
69 }
70
71
get_data_column(startree_t * s,const char * colname,const int * inds,int N,tfits_type tt)72 static void* get_data_column(startree_t* s, const char* colname, const int* inds, int N, tfits_type tt) {
73 fitstable_t* table;
74 void* arr;
75 if (N == 0) {
76 logmsg("Warning: zero stars (elements) in your request for data column \"%s\"\n", colname);
77 return NULL;
78 }
79 table = startree_get_tagalong(s);
80 if (!table) {
81 ERROR("No tag-along data found");
82 return NULL;
83 }
84 arr = fitstable_read_column_inds(table, colname, tt, inds, N);
85 if (!arr) {
86 ERROR("Failed to read tag-along data column \"%s\"", colname);
87 return NULL;
88 }
89 return arr;
90 }
91
92
startree_get_data_column(startree_t * s,const char * colname,const int * inds,int N)93 double* startree_get_data_column(startree_t* s, const char* colname, const int* inds, int N) {
94 return get_data_column(s, colname, inds, N, fitscolumn_double_type());
95 }
96
startree_get_data_column_int64(startree_t * s,const char * colname,const int * inds,int N)97 int64_t* startree_get_data_column_int64(startree_t* s, const char* colname, const int* inds, int N) {
98 return get_data_column(s, colname, inds, N, fitscolumn_i64_type());
99 }
100
startree_get_data_column_array(startree_t * s,const char * colname,const int * indices,int N,int * arraysize)101 double* startree_get_data_column_array(startree_t* s, const char* colname, const int* indices, int N, int* arraysize) {
102 fitstable_t* table;
103 tfits_type dubl = fitscolumn_double_type();
104 double* arr;
105 table = startree_get_tagalong(s);
106 if (!table) {
107 ERROR("No tag-along data found");
108 return NULL;
109 }
110 arr = fitstable_read_column_array_inds(table, colname, dubl, indices, N, arraysize);
111 if (!arr) {
112 ERROR("Failed to read tag-along data");
113 return NULL;
114 }
115 return arr;
116 }
117
startree_free_data_column(startree_t * s,double * d)118 void startree_free_data_column(startree_t* s, double* d) {
119 free(d);
120 }
121
startree_search_for_radec(const startree_t * s,double ra,double dec,double radius,double ** xyzresults,double ** radecresults,int ** starinds,int * nresults)122 void startree_search_for_radec(const startree_t* s, double ra, double dec, double radius,
123 double** xyzresults, double** radecresults,
124 int** starinds, int* nresults) {
125 double xyz[3];
126 double r2;
127 radecdeg2xyzarr(ra, dec, xyz);
128 r2 = deg2distsq(radius);
129 startree_search_for(s, xyz, r2, xyzresults, radecresults, starinds, nresults);
130 }
131
startree_search_for(const startree_t * s,const double * xyzcenter,double radius2,double ** xyzresults,double ** radecresults,int ** starinds,int * nresults)132 void startree_search_for(const startree_t* s, const double* xyzcenter, double radius2,
133 double** xyzresults, double** radecresults,
134 int** starinds, int* nresults) {
135 kdtree_qres_t* res = NULL;
136 int opts;
137 double* xyz;
138 int i, N;
139
140 opts = KD_OPTIONS_SMALL_RADIUS;
141 if (xyzresults || radecresults)
142 opts |= KD_OPTIONS_RETURN_POINTS;
143
144 res = kdtree_rangesearch_options(s->tree, xyzcenter, radius2, opts);
145
146 if (!res || !res->nres) {
147 if (xyzresults)
148 *xyzresults = NULL;
149 if (radecresults)
150 *radecresults = NULL;
151 if (starinds)
152 *starinds = NULL;
153 *nresults = 0;
154 if (res)
155 kdtree_free_query(res);
156 return;
157 }
158
159 xyz = res->results.d;
160 N = res->nres;
161 *nresults = N;
162
163 if (radecresults) {
164 *radecresults = malloc(N * 2 * sizeof(double));
165 for (i=0; i<N; i++)
166 xyzarr2radecdegarr(xyz + i*3, (*radecresults) + i*2);
167 }
168 if (xyzresults) {
169 // Steal the results array.
170 *xyzresults = xyz;
171 res->results.d = NULL;
172 }
173 if (starinds) {
174 *starinds = malloc(res->nres * sizeof(int));
175 for (i=0; i<N; i++)
176 (*starinds)[i] = res->inds[i];
177 }
178 kdtree_free_query(res);
179 }
180
181
startree_search(const startree_t * s,const double * xyzcenter,double radius2,double ** xyzresults,double ** radecresults,int * nresults)182 void startree_search(const startree_t* s, const double* xyzcenter, double radius2,
183 double** xyzresults, double** radecresults, int* nresults) {
184 startree_search_for(s, xyzcenter, radius2, xyzresults, radecresults, NULL, nresults);
185 }
186
startree_N(const startree_t * s)187 int startree_N(const startree_t* s) {
188 return s->tree->ndata;
189 }
190
startree_nodes(const startree_t * s)191 int startree_nodes(const startree_t* s) {
192 return s->tree->nnodes;
193 }
194
startree_D(const startree_t * s)195 int startree_D(const startree_t* s) {
196 return s->tree->ndim;
197 }
198
startree_header(const startree_t * s)199 qfits_header* startree_header(const startree_t* s) {
200 return s->header;
201 }
202
get_chunks(startree_t * s,il * wordsizes)203 static bl* get_chunks(startree_t* s, il* wordsizes) {
204 bl* chunks = bl_new(4, sizeof(fitsbin_chunk_t));
205 fitsbin_chunk_t chunk;
206 kdtree_t* kd = s->tree;
207
208 fitsbin_chunk_init(&chunk);
209 chunk.tablename = "sweep";
210 chunk.forced_type = fitscolumn_u8_type();
211 chunk.itemsize = sizeof(uint8_t);
212 chunk.nrows = kd->ndata;
213 chunk.data = s->sweep;
214 chunk.userdata = &(s->sweep);
215 chunk.required = FALSE;
216 bl_append(chunks, &chunk);
217 if (wordsizes)
218 il_append(wordsizes, sizeof(uint8_t));
219
220 fitsbin_chunk_clean(&chunk);
221 return chunks;
222 }
223
my_open(const char * fn,anqfits_t * fits)224 static startree_t* my_open(const char* fn, anqfits_t* fits) {
225 struct timeval tv1, tv2;
226 startree_t* s;
227 bl* chunks;
228 int i;
229 kdtree_fits_t* io;
230 char* treename = STARTREE_NAME;
231 const char* thefn = fn;
232
233 assert(fn || fits);
234
235 if (!thefn)
236 thefn = fits->filename;
237
238 s = startree_alloc();
239 if (!s)
240 return NULL;
241
242 gettimeofday(&tv1, NULL);
243 if (fn)
244 io = kdtree_fits_open(fn);
245 else
246 io = kdtree_fits_open_fits(fits);
247
248 gettimeofday(&tv2, NULL);
249 debug("kdtree_fits_open() took %g ms\n", millis_between(&tv1, &tv2));
250 if (!io) {
251 ERROR("Failed to open FITS file \"%s\"", thefn);
252 goto bailout;
253 }
254
255 gettimeofday(&tv1, NULL);
256 if (!kdtree_fits_contains_tree(io, treename))
257 treename = NULL;
258 gettimeofday(&tv2, NULL);
259 debug("kdtree_fits_contains_tree() took %g ms\n", millis_between(&tv1, &tv2));
260
261 gettimeofday(&tv1, NULL);
262 s->tree = kdtree_fits_read_tree(io, treename, &s->header);
263 gettimeofday(&tv2, NULL);
264 debug("kdtree_fits_read_tree() took %g ms\n", millis_between(&tv1, &tv2));
265 if (!s->tree) {
266 ERROR("Failed to read kdtree from file \"%s\"", thefn);
267 goto bailout;
268 }
269
270 // Check the tree dimensionality.
271 // (because code trees can be confused...)
272 if (s->tree->ndim != 3) {
273 logverb("File %s contains a kd-tree with dim %i (not 3), named %s\n",
274 thefn, s->tree->ndim, treename);
275 s->tree->io = NULL;
276 goto bailout;
277 }
278
279 gettimeofday(&tv1, NULL);
280 chunks = get_chunks(s, NULL);
281 for (i=0; i<bl_size(chunks); i++) {
282 fitsbin_chunk_t* chunk = bl_access(chunks, i);
283 void** dest = chunk->userdata;
284 kdtree_fits_read_chunk(io, chunk);
285 *dest = chunk->data;
286 }
287 bl_free(chunks);
288 gettimeofday(&tv2, NULL);
289 debug("reading chunks took %g ms\n", millis_between(&tv1, &tv2));
290
291 // kdtree_fits_t is a typedef of fitsbin_t
292 fitsbin_close_fd(io);
293
294 return s;
295
296 bailout:
297 kdtree_fits_io_close(io);
298 startree_close(s);
299 return NULL;
300 }
301
startree_open_fits(anqfits_t * fits)302 startree_t* startree_open_fits(anqfits_t* fits) {
303 return my_open(NULL, fits);
304 }
305
startree_open(const char * fn)306 startree_t* startree_open(const char* fn) {
307 return my_open(fn, NULL);
308 }
309
310 /*
311 uint64_t startree_get_starid(const startree_t* s, int ind) {
312 if (!s->starids)
313 return 0;
314 return s->starids[ind];
315 }
316 */
startree_close(startree_t * s)317 int startree_close(startree_t* s) {
318 if (!s) return 0;
319 if (s->inverse_perm)
320 free(s->inverse_perm);
321 if (s->header)
322 qfits_header_destroy(s->header);
323 if (s->tree) {
324 if (s->writing) {
325 free(s->tree->data.any);
326 s->tree->data.any = NULL;
327 kdtree_free(s->tree);
328 free(s->sweep);
329 }
330 else
331 kdtree_fits_close(s->tree);
332 }
333 if (s->tagalong)
334 fitstable_close(s->tagalong);
335 free(s);
336 return 0;
337 }
338
get_tagalong(startree_t * s,anbool report_errs)339 static fitstable_t* get_tagalong(startree_t* s, anbool report_errs) {
340 char* fn;
341 int next;
342 int i;
343 int ext = -1;
344 fitstable_t* tag;
345
346 if (!s->tree->io)
347 return NULL;
348 fn = fitsbin_get_filename(s->tree->io);
349 if (!fn) {
350 if (report_errs)
351 ERROR("No filename");
352 return NULL;
353 }
354 tag = fitstable_open(fn);
355 if (!tag) {
356 if (report_errs)
357 ERROR("Failed to open FITS table from %s", fn);
358 return NULL;
359 }
360 next = fitstable_n_extensions(tag);
361 for (i=1; i<next; i++) {
362 char* type;
363 anbool eq;
364 const qfits_header* hdr;
365 hdr = anqfits_get_header_const(tag->anq, i);
366 if (!hdr) {
367 if (report_errs)
368 ERROR("Failed to read FITS header for ext %i in %s", i, fn);
369 return NULL;
370 }
371 type = fits_get_dupstring(hdr, "AN_FILE");
372 eq = streq(type, AN_FILETYPE_TAGALONG);
373 free(type);
374 if (!eq)
375 continue;
376 ext = i;
377 break;
378 }
379 if (ext == -1) {
380 if (report_errs)
381 ERROR("Failed to find a FITS header with the card AN_FILE = TAGALONG");
382 return NULL;
383 }
384 fitstable_open_extension(tag, ext);
385 return tag;
386 }
387
startree_get_tagalong(startree_t * s)388 fitstable_t* startree_get_tagalong(startree_t* s) {
389 if (s->tagalong)
390 return s->tagalong;
391 s->tagalong = get_tagalong(s, TRUE);
392 return s->tagalong;
393 }
394
startree_has_tagalong(startree_t * s)395 anbool startree_has_tagalong(startree_t* s) {
396 return (startree_get_tagalong(s) != NULL);
397 }
398
Ndata(const startree_t * s)399 static int Ndata(const startree_t* s) {
400 return s->tree->ndata;
401 }
402
startree_check_inverse_perm(startree_t * s)403 int startree_check_inverse_perm(startree_t* s) {
404 // ensure that each value appears exactly once.
405 int i, N;
406 uint8_t* counts;
407 N = Ndata(s);
408 counts = calloc(Ndata(s), sizeof(uint8_t));
409 for (i=0; i<N; i++) {
410 assert(s->inverse_perm[i] >= 0);
411 assert(s->inverse_perm[i] < N);
412 counts[s->inverse_perm[i]]++;
413 }
414 for (i=0; i<N; i++) {
415 assert(counts[i] == 1);
416 }
417 free(counts); //# Modified by Robert Lancaster for the StellarSolver Internal Library, to prevent leak
418 return 0;
419 }
420
startree_compute_inverse_perm(startree_t * s)421 void startree_compute_inverse_perm(startree_t* s) {
422 if (s->inverse_perm)
423 return;
424 // compute inverse permutation vector.
425 s->inverse_perm = malloc(Ndata(s) * sizeof(int));
426 if (!s->inverse_perm) {
427 fprintf(stderr, "Failed to allocate star kdtree inverse permutation vector.\n");
428 return;
429 }
430 #ifndef NDEBUG
431 {
432 int i;
433 for (i=0; i<Ndata(s); i++)
434 s->inverse_perm[i] = -1;
435 }
436 #endif
437 kdtree_inverse_permutation(s->tree, s->inverse_perm);
438 #ifndef NDEBUG
439 {
440 int i;
441 for (i=0; i<Ndata(s); i++)
442 assert(s->inverse_perm[i] != -1);
443 }
444 #endif
445 }
446
startree_get_cut_nside(const startree_t * s)447 int startree_get_cut_nside(const startree_t* s) {
448 return qfits_header_getint(s->header, "CUTNSIDE", -1);
449 }
450
startree_get_cut_nsweeps(const startree_t * s)451 int startree_get_cut_nsweeps(const startree_t* s) {
452 return qfits_header_getint(s->header, "CUTNSWEP", -1);
453 }
454
startree_get_cut_dedup(const startree_t * s)455 double startree_get_cut_dedup(const startree_t* s) {
456 return qfits_header_getdouble(s->header, "CUTDEDUP", 0.0);
457 }
458
startree_get_cut_band(const startree_t * s)459 char* startree_get_cut_band(const startree_t* s) {
460 static char* bands[] = { "R", "B", "J" };
461 int i;
462 char* str = fits_get_dupstring(s->header, "CUTBAND");
463 char* rtn = NULL;
464 if (!str)
465 return NULL;
466 for (i=0; i<sizeof(bands) / sizeof(char*); i++) {
467 if (streq(str, bands[i])) {
468 rtn = bands[i];
469 break;
470 }
471 }
472 free(str);
473 return rtn;
474 }
475
startree_get_cut_margin(const startree_t * s)476 int startree_get_cut_margin(const startree_t* s) {
477 return qfits_header_getint(s->header, "CUTMARG", -1);
478 }
479
startree_get_jitter(const startree_t * s)480 double startree_get_jitter(const startree_t* s) {
481 return qfits_header_getdouble(s->header, "JITTER", 0.0);
482 }
483
startree_set_jitter(startree_t * s,double jitter_arcsec)484 void startree_set_jitter(startree_t* s, double jitter_arcsec) {
485 fits_header_set_double(s->header, "JITTER", jitter_arcsec, "Positional error of stars [arcsec]");
486 }
487
startree_get_sweep(const startree_t * s,int ind)488 int startree_get_sweep(const startree_t* s, int ind) {
489 if (ind < 0 || ind >= Ndata(s) || !s->sweep)
490 return -1;
491 return s->sweep[ind];
492 }
493
startree_get(startree_t * s,int starid,double * posn)494 int startree_get(startree_t* s, int starid, double* posn) {
495 if (s->tree->perm && !s->inverse_perm) {
496 startree_compute_inverse_perm(s);
497 if (!s->inverse_perm)
498 return -1;
499 }
500 if (starid >= Ndata(s)) {
501 fprintf(stderr, "Invalid star ID: %u >= %u.\n", starid, Ndata(s));
502 assert(0);
503 return -1;
504 }
505 if (s->inverse_perm) {
506 kdtree_copy_data_double(s->tree, s->inverse_perm[starid], 1, posn);
507 } else {
508 kdtree_copy_data_double(s->tree, starid, 1, posn);
509 }
510 return 0;
511 }
512
startree_get_radec(startree_t * s,int starid,double * ra,double * dec)513 int startree_get_radec(startree_t* s, int starid, double* ra, double* dec) {
514 double xyz[3];
515 int rtn;
516 rtn = startree_get(s, starid, xyz);
517 if (rtn)
518 return rtn;
519 xyzarr2radecdeg(xyz, ra, dec);
520 return rtn;
521 }
522
startree_new()523 startree_t* startree_new() {
524 startree_t* s = startree_alloc();
525 s->header = qfits_header_default();
526 if (!s->header) {
527 fprintf(stderr, "Failed to create a qfits header for star kdtree.\n");
528 free(s);
529 return NULL;
530 }
531 qfits_header_add(s->header, "AN_FILE", AN_FILETYPE_STARTREE, "This file is a star kdtree.", NULL);
532 s->writing = TRUE;
533 return s;
534 }
535
write_to_file(startree_t * s,const char * fn,anbool flipped,FILE * fid)536 static int write_to_file(startree_t* s, const char* fn, anbool flipped,
537 FILE* fid) {
538 bl* chunks;
539 il* wordsizes = NULL;
540 int i;
541 kdtree_fits_t* io = NULL;
542
543 // just haven't bothered...
544 assert(!(flipped && fid));
545
546 int status = 0; //# Modified by Robert Lancaster for the StellarSolver Internal Library, to prevent leak
547
548 if (fn) {
549 io = kdtree_fits_open_for_writing(fn);
550 if (!io) {
551 ERROR("Failed to open file \"%s\" for writing kdtree", fn);
552 status = -1;
553 goto exit; //# Modified by Robert Lancaster for the StellarSolver Internal Library, to prevent leak
554 }
555 }
556 if (flipped) {
557 if (kdtree_fits_write_tree_flipped(io, s->tree, s->header)) {
558 ERROR("Failed to write (flipped) kdtree to file \"%s\"", fn);
559 status = -1;
560 goto exit; //# Modified by Robert Lancaster for the StellarSolver Internal Library, to prevent leak
561 }
562 } else {
563 if (fid) {
564 if (kdtree_fits_append_tree_to(s->tree, s->header, fid)) {
565 ERROR("Failed to write star kdtree");
566 status = -1;
567 goto exit; //# Modified by Robert Lancaster for the StellarSolver Internal Library, to prevent leak
568 }
569 } else {
570 if (kdtree_fits_write_tree(io, s->tree, s->header)) {
571 ERROR("Failed to write kdtree to file \"%s\"", fn);
572 status = -1;
573 goto exit; //# Modified by Robert Lancaster for the StellarSolver Internal Library, to prevent leak
574 }
575 }
576 }
577
578 if (flipped)
579 wordsizes = il_new(4);
580
581 chunks = get_chunks(s, wordsizes);
582 for (i=0; i<bl_size(chunks); i++) {
583 fitsbin_chunk_t* chunk = bl_access(chunks, i);
584 if (!chunk->data)
585 continue;
586 if (flipped)
587 kdtree_fits_write_chunk_flipped(io, chunk, il_get(wordsizes, i));
588 else {
589 if (fid) {
590 kdtree_fits_write_chunk_to(chunk, fid);
591 } else {
592 kdtree_fits_write_chunk(io, chunk);
593 }
594 }
595 fitsbin_chunk_clean(chunk);
596 }
597 bl_free(chunks);
598
599 if (flipped)
600 il_free(wordsizes);
601
602 exit: //# Modified by Robert Lancaster for the StellarSolver Internal Library, to prevent leak
603 if (io)
604 kdtree_fits_io_close(io);
605 return 0;
606 }
607
startree_write_to_file(startree_t * s,const char * fn)608 int startree_write_to_file(startree_t* s, const char* fn) {
609 return write_to_file(s, fn, FALSE, NULL);
610 }
611
startree_write_to_file_flipped(startree_t * s,const char * fn)612 int startree_write_to_file_flipped(startree_t* s, const char* fn) {
613 return write_to_file(s, fn, TRUE, NULL);
614 }
615
startree_append_to(startree_t * s,FILE * fid)616 int startree_append_to(startree_t* s, FILE* fid) {
617 return write_to_file(s, NULL, FALSE, fid);
618 }
619
620