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 <unistd.h>
7 #include <errno.h>
8 #include <string.h>
9 #include <stdio.h>
10 #include <math.h>
11 #include <assert.h>
12
13 #include "os-features.h"
14 #include "catalog.h"
15 #include "fitsioutils.h"
16 #include "ioutils.h"
17 #include "starutil.h"
18 #include "mathutil.h"
19 #include "errors.h"
20
21 #define CHUNK_XYZ 0
22 #define CHUNK_MAG 1
23 #define CHUNK_MAG_ERR 2
24 #define CHUNK_SIG 3
25 #define CHUNK_PM 4
26 #define CHUNK_SIGPM 5
27 #define CHUNK_STARID 6
28
xyz_chunk(catalog * cat)29 static fitsbin_chunk_t* xyz_chunk(catalog* cat) {
30 return fitsbin_get_chunk(cat->fb, CHUNK_XYZ);
31 }
mag_chunk(catalog * cat)32 static fitsbin_chunk_t* mag_chunk(catalog* cat) {
33 return fitsbin_get_chunk(cat->fb, CHUNK_MAG);
34 }
mag_err_chunk(catalog * cat)35 static fitsbin_chunk_t* mag_err_chunk(catalog* cat) {
36 return fitsbin_get_chunk(cat->fb, CHUNK_MAG_ERR);
37 }
get_chunk(catalog * cat,int i)38 static fitsbin_chunk_t* get_chunk(catalog* cat, int i) {
39 return fitsbin_get_chunk(cat->fb, i);
40 }
41
callback_read_header(fitsbin_t * fb,fitsbin_chunk_t * chunk)42 static int callback_read_header(fitsbin_t* fb, fitsbin_chunk_t* chunk) {
43 qfits_header* primheader = fitsbin_get_primary_header(fb);
44 catalog* cat = chunk->userdata;
45 cat->numstars = qfits_header_getint(primheader, "NSTARS", -1);
46 cat->healpix = qfits_header_getint(primheader, "HEALPIX", -1);
47 cat->hpnside = qfits_header_getint(primheader, "HPNSIDE", 1);
48 if (fits_check_endian(primheader)) {
49 ERROR("Catalog file was written with wrong endianness");
50 return -1;
51 }
52 if (cat->numstars == -1) {
53 ERROR("Couldn't find NSTARS header in catalog file.");
54 return -1;
55 }
56
57 chunk->nrows = cat->numstars;
58 return 0;
59 }
60
callback_read_tagalong(fitsbin_t * fb,fitsbin_chunk_t * chunk)61 static int callback_read_tagalong(fitsbin_t* fb, fitsbin_chunk_t* chunk) {
62 catalog* cat = chunk->userdata;
63 chunk->nrows = cat->numstars;
64 return 0;
65 }
66
new_catalog(const char * fn,anbool writing)67 static catalog* new_catalog(const char* fn, anbool writing) {
68 catalog* cat;
69 fitsbin_chunk_t chunk;
70
71 cat = calloc(1, sizeof(catalog));
72 if (!cat) {
73 fprintf(stderr, "catalog_open: malloc failed.\n");
74 }
75
76 if (writing)
77 cat->fb = fitsbin_open_for_writing(fn);
78 else
79 cat->fb = fitsbin_open(fn);
80 if (!cat->fb) {
81 ERROR("Failed to create fitsbin");
82 return NULL;
83 }
84
85 memset(&chunk, 0, sizeof(fitsbin_chunk_t));
86
87 // NOTE -- the order these are added MUST match the CHUNK_XYZ, CHUNK_MAG, etc
88 // ordering.
89
90 // Star positions
91 fitsbin_chunk_init(&chunk);
92 chunk.tablename = "xyz";
93 chunk.required = 1;
94 chunk.callback_read_header = callback_read_header;
95 chunk.userdata = cat;
96 chunk.itemsize = DIM_STARS * sizeof(double);
97 fitsbin_add_chunk(cat->fb, &chunk);
98 fitsbin_chunk_reset(&chunk);
99
100 // Star magnitudes
101 chunk.tablename = "mag";
102 chunk.required = 0;
103 chunk.callback_read_header = callback_read_tagalong;
104 chunk.userdata = cat;
105 chunk.itemsize = sizeof(float);
106 fitsbin_add_chunk(cat->fb, &chunk);
107 fitsbin_chunk_reset(&chunk);
108
109 // Star magnitude errors
110 chunk.tablename = "mag_err";
111 chunk.required = 0;
112 chunk.callback_read_header = callback_read_tagalong;
113 chunk.userdata = cat;
114 chunk.itemsize = sizeof(float);
115 fitsbin_add_chunk(cat->fb, &chunk);
116 fitsbin_chunk_reset(&chunk);
117
118 // Sigmas
119 chunk.tablename = "sigma_radec";
120 chunk.required = 0;
121 chunk.callback_read_header = callback_read_tagalong;
122 chunk.userdata = cat;
123 chunk.itemsize = 2 * sizeof(float);
124 fitsbin_add_chunk(cat->fb, &chunk);
125 fitsbin_chunk_reset(&chunk);
126
127 // Motions
128 chunk.tablename = "proper_motion";
129 chunk.required = 0;
130 chunk.callback_read_header = callback_read_tagalong;
131 chunk.userdata = cat;
132 chunk.itemsize = 2 * sizeof(float);
133 fitsbin_add_chunk(cat->fb, &chunk);
134 fitsbin_chunk_reset(&chunk);
135
136 // Sigma Motions
137 chunk.tablename = "sigma_pm";
138 chunk.required = 0;
139 chunk.callback_read_header = callback_read_tagalong;
140 chunk.userdata = cat;
141 chunk.itemsize = 2 * sizeof(float);
142 fitsbin_add_chunk(cat->fb, &chunk);
143 fitsbin_chunk_reset(&chunk);
144
145 // Ids
146 chunk.tablename = "starid";
147 chunk.required = 0;
148 chunk.callback_read_header = callback_read_tagalong;
149 chunk.userdata = cat;
150 chunk.itemsize = sizeof(uint64_t);
151 fitsbin_add_chunk(cat->fb, &chunk);
152 fitsbin_chunk_clean(&chunk);
153
154 return cat;
155 }
156
catalog_has_mag(const catalog * cat)157 anbool catalog_has_mag(const catalog* cat) {
158 return (cat->mag != NULL);
159 }
160
catalog_open(char * fn)161 catalog* catalog_open(char* fn) {
162 catalog* cat = NULL;
163
164 cat = new_catalog(fn, FALSE);
165 if (!cat)
166 goto bailout;
167 if (fitsbin_read(cat->fb)) {
168 ERROR("catalog: fitsbin_read() failed");
169 goto bailout;
170 }
171 cat->stars = xyz_chunk(cat)->data;
172 cat->mag = mag_chunk(cat)->data;
173 cat->mag_err = mag_err_chunk(cat)->data;
174 cat->sigma_radec = get_chunk(cat, CHUNK_SIG )->data;
175 cat->proper_motion = get_chunk(cat, CHUNK_PM )->data;
176 cat->sigma_pm = get_chunk(cat, CHUNK_SIGPM )->data;
177 cat->starids = get_chunk(cat, CHUNK_STARID)->data;
178 return cat;
179
180 bailout:
181 if (cat)
182 catalog_close(cat);
183 return NULL;
184 }
185
add_default_header_vals(catalog * cat)186 static void add_default_header_vals(catalog* cat) {
187 qfits_header* hdr;
188 hdr = catalog_get_header(cat);
189 qfits_header_add(hdr, "AN_FILE", "OBJS", "This file has a list of object positions.", NULL);
190 fits_add_endian(hdr);
191 fits_add_double_size(hdr);
192 qfits_header_add(hdr, "NSTARS", "0", "Number of stars in this file.", NULL);
193 qfits_header_add(hdr, "HEALPIX", "-1", "Healpix covered by this catalog, with Nside=HPNSIDE", NULL);
194 qfits_header_add(hdr, "HPNSIDE", "-1", "Nside of HEALPIX.", NULL);
195 qfits_header_add(hdr, "COMMENT", "This is a flat array of XYZ for each catalog star.", NULL, NULL);
196 qfits_header_add(hdr, "COMMENT", " (ie, star position on the unit sphere)", NULL, NULL);
197 qfits_header_add(hdr, "COMMENT", " (stored as three native-{endian,size} doubles)", NULL, NULL);
198 }
199
catalog_open_for_writing(char * fn)200 catalog* catalog_open_for_writing(char* fn) {
201 catalog* cat;
202
203 cat = new_catalog(fn, TRUE);
204 if (!cat)
205 goto bailout;
206 cat->hpnside = 1;
207 add_default_header_vals(cat);
208 return cat;
209
210 bailout:
211 if (cat)
212 catalog_close(cat);
213 return NULL;
214 }
215
catalog_get_header(catalog * cat)216 qfits_header* catalog_get_header(catalog* cat) {
217 return fitsbin_get_primary_header(cat->fb);
218 }
219
catalog_write_header(catalog * cat)220 int catalog_write_header(catalog* cat) {
221 fitsbin_t* fb = cat->fb;
222
223 if (fitsbin_write_primary_header(fb) ||
224 fitsbin_write_chunk_header(fb, xyz_chunk(cat))) {
225 ERROR("Failed to write catalog header");
226 return -1;
227 }
228 return 0;
229 }
230
catalog_fix_header(catalog * cat)231 int catalog_fix_header(catalog* cat) {
232 qfits_header* hdr;
233 fitsbin_t* fb = cat->fb;
234
235 hdr = catalog_get_header(cat);
236 // fill in the real values...
237 fits_header_mod_int(hdr, "NSTARS", cat->numstars, "Number of stars.");
238 fits_header_mod_int(hdr, "HEALPIX", cat->healpix, "Healpix covered by this catalog, with Nside=HPNSIDE");
239 fits_header_mod_int(hdr, "HPNSIDE", cat->hpnside, "Nside of HEALPIX.");
240
241 if (fitsbin_fix_primary_header(fb) ||
242 fitsbin_fix_chunk_header(fb, xyz_chunk(cat))) {
243 ERROR("Failed to fix catalog header");
244 return -1;
245 }
246 return 0;
247 }
248
catalog_get_base(catalog * cat)249 double* catalog_get_base(catalog* cat) {
250 return cat->stars;
251 }
252
catalog_get_star(catalog * cat,int sid)253 double* catalog_get_star(catalog* cat, int sid) {
254 if (sid >= cat->numstars) {
255 fflush(stdout);
256 fprintf(stderr, "catalog: asked for star %u, but catalog size is only %u.\n",
257 sid, cat->numstars);
258 return NULL;
259 }
260 return cat->stars + sid * 3;
261 }
262
catalog_write_star(catalog * cat,double * star)263 int catalog_write_star(catalog* cat, double* star) {
264 if (fitsbin_write_item(cat->fb, xyz_chunk(cat), star)) {
265 fprintf(stderr, "Failed to write catalog data to file %s: %s\n",
266 cat->fb->filename, strerror(errno));
267 return -1;
268 }
269 cat->numstars++;
270 return 0;
271 }
272
write_floats(catalog * cat,int chunknum,const char * name,fl * list,int nperstar)273 int write_floats(catalog* cat, int chunknum,
274 const char* name, fl* list, int nperstar) {
275 int i;
276 int B = 1000;
277 fitsbin_chunk_t* chunk = get_chunk(cat, chunknum);
278 if (!list || (fl_size(list) != cat->numstars * nperstar)) {
279 ERROR("Number of %ss (%i) doesn't match number of stars (%i)",
280 name, list ? fl_size(list) : 0, cat->numstars);
281 return -1;
282 }
283
284 if (fitsbin_write_chunk_header(cat->fb, chunk)) {
285 ERROR("Failed to write %ss header", name);
286 return -1;
287 }
288 for (i=0; i<cat->numstars; i+=B) {
289 float data[nperstar * B];
290 int n = MIN(i+B, cat->numstars) - i;
291 fl_copy(list, i*nperstar, nperstar*n, data);
292 if (fitsbin_write_items(cat->fb, chunk, data, n)) {
293 ERROR("Failed to write %s for stars %i to %i", name, i, i+n-1);
294 return -1;
295 }
296 }
297 if (fitsbin_fix_chunk_header(cat->fb, chunk)) {
298 ERROR("Failed to fix %ss header", name);
299 return -1;
300 }
301 return 0;
302 }
303
catalog_write_mags(catalog * cat)304 int catalog_write_mags(catalog* cat) {
305 return write_floats(cat, CHUNK_MAG, "magnitude", cat->maglist, 1);
306 }
307
catalog_write_mag_errs(catalog * cat)308 int catalog_write_mag_errs(catalog* cat) {
309 return write_floats(cat, CHUNK_MAG_ERR, "magnitude errors", cat->magerrlist, 1);
310 }
311
catalog_write_sigmas(catalog * cat)312 int catalog_write_sigmas(catalog* cat) {
313 return write_floats(cat, CHUNK_SIG, "sigma", cat->siglist, 2);
314 }
315
catalog_write_pms(catalog * cat)316 int catalog_write_pms(catalog* cat) {
317 return write_floats(cat, CHUNK_PM, "proper motion", cat->pmlist, 2);
318 }
319
catalog_write_sigma_pms(catalog * cat)320 int catalog_write_sigma_pms(catalog* cat) {
321 return write_floats(cat, CHUNK_SIGPM, "sigma proper motion", cat->sigpmlist, 2);
322 }
323
catalog_write_ids(catalog * cat)324 int catalog_write_ids(catalog* cat) {
325 int i;
326 char* name = "id";
327 int chunknum = CHUNK_STARID;
328 fitsbin_chunk_t* chunk = get_chunk(cat, chunknum);
329 if (!cat->idlist || (bl_size(cat->idlist) != cat->numstars)) {
330 ERROR("Number of %ss (%i) doesn't match number of stars (%i)",
331 name, cat->idlist ? bl_size(cat->idlist) : 0, cat->numstars);
332 return -1;
333 }
334
335 if (fitsbin_write_chunk_header(cat->fb, chunk)) {
336 ERROR("Failed to write %ss header", name);
337 return -1;
338 }
339 for (i=0; i<cat->numstars; i++)
340 if (fitsbin_write_item(cat->fb, chunk, bl_access(cat->idlist, i))) {
341 ERROR("Failed to write %s for star %i", name, i);
342 return -1;
343 }
344 if (fitsbin_fix_chunk_header(cat->fb, chunk)) {
345 ERROR("Failed to fix %ss header", name);
346 return -1;
347 }
348 return 0;
349 }
350
addfloat(fl ** list,float val)351 void addfloat(fl** list, float val) {
352 if (!(*list))
353 *list = fl_new(256);
354 fl_append(*list, val);
355 }
356
catalog_add_mag(catalog * cat,float mag)357 void catalog_add_mag(catalog* cat, float mag) {
358 addfloat(&(cat->maglist), mag);
359 }
catalog_add_mag_err(catalog * cat,float magerr)360 void catalog_add_mag_err(catalog* cat, float magerr) {
361 addfloat(&(cat->magerrlist), magerr);
362 }
catalog_add_sigmas(catalog * cat,float sra,float sdec)363 void catalog_add_sigmas(catalog* cat, float sra, float sdec) {
364 addfloat(&(cat->siglist), sra);
365 addfloat(&(cat->siglist), sdec);
366 }
catalog_add_pms(catalog * cat,float sra,float sdec)367 void catalog_add_pms(catalog* cat, float sra, float sdec) {
368 addfloat(&(cat->pmlist), sra);
369 addfloat(&(cat->pmlist), sdec);
370 }
catalog_add_sigma_pms(catalog * cat,float sra,float sdec)371 void catalog_add_sigma_pms(catalog* cat, float sra, float sdec) {
372 addfloat(&(cat->sigpmlist), sra);
373 addfloat(&(cat->sigpmlist), sdec);
374 }
catalog_add_id(catalog * cat,uint64_t id)375 void catalog_add_id(catalog* cat, uint64_t id) {
376 if (!cat->idlist)
377 cat->idlist = bl_new(256, sizeof(uint64_t));
378 bl_append(cat->idlist, &id);
379 }
380
catalog_close(catalog * cat)381 int catalog_close(catalog* cat) {
382 int rtn;
383 if (!cat) return 0;
384 rtn = fitsbin_close(cat->fb);
385 fl_free(cat->maglist);
386 fl_free(cat->magerrlist);
387 fl_free(cat->siglist);
388 fl_free(cat->pmlist);
389 fl_free(cat->sigpmlist);
390 bl_free(cat->idlist);
391 free(cat);
392 return rtn;
393 }
394
395