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 <math.h>
8 #include <string.h>
9 #include <assert.h>
10
11 #include "quadfile.h"
12 #include "qfits_header.h"
13 #include "fitsioutils.h"
14 #include "starutil.h"
15 #include "ioutils.h"
16 #include "errors.h"
17 #include "an-endian.h"
18
19 #define CHUNK_QUADS 0
20
quads_chunk(quadfile_t * qf)21 static fitsbin_chunk_t* quads_chunk(quadfile_t* qf) {
22 return fitsbin_get_chunk(qf->fb, CHUNK_QUADS);
23 }
24
callback_read_header(fitsbin_t * fb,fitsbin_chunk_t * chunk)25 static int callback_read_header(fitsbin_t* fb, fitsbin_chunk_t* chunk) {
26 qfits_header* primheader = fitsbin_get_primary_header(fb);
27 quadfile_t* qf = chunk->userdata;
28
29 qf->dimquads = qfits_header_getint(primheader, "DIMQUADS", 4);
30 qf->numquads = qfits_header_getint(primheader, "NQUADS", -1);
31 qf->numstars = qfits_header_getint(primheader, "NSTARS", -1);
32 qf->index_scale_upper = qfits_header_getdouble(primheader, "SCALE_U", -1.0);
33 qf->index_scale_lower = qfits_header_getdouble(primheader, "SCALE_L", -1.0);
34 qf->indexid = qfits_header_getint(primheader, "INDEXID", 0);
35 qf->healpix = qfits_header_getint(primheader, "HEALPIX", -1);
36 qf->hpnside = qfits_header_getint(primheader, "HPNSIDE", 1);
37
38 if ((qf->numquads == -1) || (qf->numstars == -1) ||
39 (qf->index_scale_upper == -1.0) || (qf->index_scale_lower == -1.0)) {
40 ERROR("Couldn't find NQUADS or NSTARS or SCALE_U or SCALE_L entries in FITS header");
41 return -1;
42 }
43 if (fits_check_endian(primheader)) {
44 ERROR("Quad file was written with the wrong endianness");
45 return -1;
46 }
47
48 chunk->itemsize = qf->dimquads * sizeof(uint32_t);
49 chunk->nrows = qf->numquads;
50 return 0;
51 }
52
new_quadfile(const char * fn,anqfits_t * fits,anbool writing)53 static quadfile_t* new_quadfile(const char* fn, anqfits_t* fits, anbool writing) {
54 quadfile_t* qf;
55 fitsbin_chunk_t chunk;
56 qf = calloc(1, sizeof(quadfile_t));
57 if (!qf) {
58 SYSERROR("Couldn't malloc a quadfile struct");
59 return NULL;
60 }
61 qf->healpix = -1;
62 qf->hpnside = 1;
63
64 if (writing)
65 if (fn) {
66 qf->fb = fitsbin_open_for_writing(fn);
67 } else {
68 qf->fb = fitsbin_open_in_memory();
69 }
70 else {
71 if (fits)
72 qf->fb = fitsbin_open_fits(fits);
73 else
74 qf->fb = fitsbin_open(fn);
75 }
76 if (!qf->fb) {
77 ERROR("Failed to create fitsbin");
78 return NULL;
79 }
80
81 fitsbin_chunk_init(&chunk);
82 chunk.tablename = "quads";
83 chunk.required = 1;
84 chunk.callback_read_header = callback_read_header;
85 chunk.userdata = qf;
86 fitsbin_add_chunk(qf->fb, &chunk);
87 fitsbin_chunk_clean(&chunk);
88
89 return qf;
90 }
91
quadfile_check(const quadfile_t * qf)92 int quadfile_check(const quadfile_t* qf) {
93 int q, i;
94 if (qf->dimquads < 3 || qf->dimquads > DQMAX) {
95 ERROR("Dimquads has illegal value %i", qf->dimquads);
96 return -1;
97 }
98 for (q=0; q<qf->numquads; q++) {
99 unsigned int stars[DQMAX];
100 if (quadfile_get_stars(qf, q, stars)) {
101 ERROR("Failed to get quad %i of %i", q, qf->numquads);
102 return -1;
103 }
104 for (i=0; i<qf->dimquads; i++) {
105 if (stars[i] >= qf->numstars) {
106 ERROR("Star ID %i is out of bounds: num stars %i", stars[i], qf->numstars);
107 return -1;
108 }
109 }
110 }
111 return 0;
112 }
113
quadfile_dimquads(const quadfile_t * qf)114 int quadfile_dimquads(const quadfile_t* qf) {
115 return qf->dimquads;
116 }
117
quadfile_nquads(const quadfile_t * qf)118 int quadfile_nquads(const quadfile_t* qf) {
119 return qf->numquads;
120 }
121
quadfile_get_header(const quadfile_t * qf)122 qfits_header* quadfile_get_header(const quadfile_t* qf) {
123 return fitsbin_get_primary_header(qf->fb);
124 }
125
my_open(const char * fn,anqfits_t * fits)126 static quadfile_t* my_open(const char* fn, anqfits_t* fits) {
127 quadfile_t* qf = NULL;
128 fitsbin_chunk_t* chunk;
129
130 qf = new_quadfile(fn, fits, FALSE);
131 if (!qf)
132 goto bailout;
133 if (fitsbin_read(qf->fb)) {
134 ERROR("Failed to open quads file");
135 goto bailout;
136 }
137 chunk = quads_chunk(qf);
138 qf->quadarray = chunk->data;
139
140 // close fd.
141 if (qf->fb->fid) {
142 if (fclose(qf->fb->fid)) {
143 ERROR("Failed to close quadfile FID");
144 goto bailout;
145 }
146 qf->fb->fid = NULL;
147 }
148
149 return qf;
150
151 bailout:
152 if (qf)
153 quadfile_close(qf);
154 return NULL;
155 }
156
quadfile_get_filename(const quadfile_t * qf)157 char* quadfile_get_filename(const quadfile_t* qf) {
158 return fitsbin_get_filename(qf->fb);
159 }
160
quadfile_open_fits(anqfits_t * fits)161 quadfile_t* quadfile_open_fits(anqfits_t* fits) {
162 return my_open(NULL, fits);
163 }
164
quadfile_open(const char * fn)165 quadfile_t* quadfile_open(const char* fn) {
166 return my_open(fn, NULL);
167 }
168
quadfile_close(quadfile_t * qf)169 int quadfile_close(quadfile_t* qf) {
170 int rtn;
171 if (!qf) return 0;
172 rtn = fitsbin_close(qf->fb);
173 free(qf);
174 return rtn;
175 }
176
open_for_writing(const char * fn)177 static quadfile_t* open_for_writing(const char* fn) {
178 quadfile_t* qf;
179 qfits_header* hdr;
180 qf = new_quadfile(fn, NULL, TRUE);
181 if (!qf)
182 goto bailout;
183 qf->dimquads = 4;
184 // add default values to header
185 hdr = fitsbin_get_primary_header(qf->fb);
186 fits_add_endian(hdr);
187 qfits_header_add(hdr, "AN_FILE", "QUAD", "This file lists, for each quad, its stars.", NULL);
188 qfits_header_add(hdr, "DIMQUADS", "0", "", NULL);
189 qfits_header_add(hdr, "NQUADS", "0", "", NULL);
190 qfits_header_add(hdr, "NSTARS", "0", "", NULL);
191 qfits_header_add(hdr, "SCALE_U", "0.0", "", NULL);
192 qfits_header_add(hdr, "SCALE_L", "0.0", "", NULL);
193 qfits_header_add(hdr, "INDEXID", "0", "", NULL);
194 qfits_header_add(hdr, "HEALPIX", "-1", "", NULL);
195 qfits_header_add(hdr, "HPNSIDE", "1", "", NULL);
196 fits_add_long_comment(hdr, "The first extension contains the quads "
197 "stored as %i 32-bit native-endian unsigned ints.", qf->dimquads);
198 return qf;
199
200 bailout:
201 if (qf)
202 quadfile_close(qf);
203 return NULL;
204 }
205
quadfile_open_for_writing(const char * fn)206 quadfile_t* quadfile_open_for_writing(const char* fn) {
207 if (!fn) {
208 ERROR("Non-NULL filename required");
209 return NULL;
210 }
211 return open_for_writing(fn);
212 }
213
quadfile_open_in_memory()214 quadfile_t* quadfile_open_in_memory() {
215 return open_for_writing(NULL);
216 }
217
quadfile_switch_to_reading(quadfile_t * qf)218 int quadfile_switch_to_reading(quadfile_t* qf) {
219 if (quadfile_fix_header(qf)) {
220 ERROR("Failed to fix quads header");
221 return -1;
222 }
223 if (fitsbin_switch_to_reading(qf->fb)) {
224 ERROR("Failed to switch to read mode");
225 return -1;
226 }
227 if (fitsbin_read(qf->fb)) {
228 ERROR("Failed to open quads file");
229 return -1;
230 }
231 qf->quadarray = quads_chunk(qf)->data;
232 return 0;
233 }
234
add_to_header(qfits_header * hdr,quadfile_t * qf)235 static void add_to_header(qfits_header* hdr, quadfile_t* qf) {
236 fits_header_mod_int(hdr, "DIMQUADS", qf->dimquads, "Number of stars in a quad.");
237 fits_header_mod_int(hdr, "NQUADS", qf->numquads, "Number of quads.");
238 fits_header_mod_int(hdr, "NSTARS", qf->numstars, "Number of stars.");
239 fits_header_mod_double(hdr, "SCALE_U", qf->index_scale_upper, "Upper-bound index scale (radians).");
240 fits_header_mod_double(hdr, "SCALE_L", qf->index_scale_lower, "Lower-bound index scale (radians).");
241 fits_header_mod_int(hdr, "INDEXID", qf->indexid, "Index unique ID.");
242 fits_header_mod_int(hdr, "HEALPIX", qf->healpix, "Healpix of this index.");
243 fits_header_mod_int(hdr, "HPNSIDE", qf->hpnside, "Nside of the healpixelization");
244 }
245
quadfile_write_header(quadfile_t * qf)246 int quadfile_write_header(quadfile_t* qf) {
247 fitsbin_t* fb = qf->fb;
248 fitsbin_chunk_t* chunk = quads_chunk(qf);
249 qfits_header* hdr;
250 chunk->itemsize = qf->dimquads * sizeof(uint32_t);
251 chunk->nrows = qf->numquads;
252
253 hdr = fitsbin_get_primary_header(fb);
254 add_to_header(hdr, qf);
255
256 if (fitsbin_write_primary_header(fb) ||
257 fitsbin_write_chunk_header(fb, chunk)) {
258 ERROR("Failed to write quadfile header");
259 return -1;
260 }
261 return 0;
262 }
263
quadfile_write_header_to(quadfile_t * qf,FILE * fid)264 int quadfile_write_header_to(quadfile_t* qf, FILE* fid) {
265 fitsbin_t* fb = qf->fb;
266 fitsbin_chunk_t* chunk = quads_chunk(qf);
267 qfits_header* hdr;
268 chunk->itemsize = qf->dimquads * sizeof(uint32_t);
269 chunk->nrows = qf->numquads;
270 hdr = fitsbin_get_primary_header(fb);
271 add_to_header(hdr, qf);
272
273 if (fitsbin_write_primary_header_to(fb, fid) ||
274 fitsbin_write_chunk_header_to(fb, chunk, fid)) {
275 ERROR("Failed to write quadfile header");
276 return -1;
277 }
278 return 0;
279 }
280
quadfile_write_quad(quadfile_t * qf,unsigned int * stars)281 int quadfile_write_quad(quadfile_t* qf, unsigned int* stars) {
282 uint32_t* data;
283 uint32_t ustars[qf->dimquads];
284 int i;
285 fitsbin_chunk_t* chunk = quads_chunk(qf);
286
287 if (sizeof(uint32_t) == sizeof(uint)) {
288 data = stars;
289 } else {
290 data = ustars;
291 for (i=0; i<qf->dimquads; i++)
292 ustars[i] = stars[i];
293 }
294 if (fitsbin_write_item(qf->fb, chunk, data)) {
295 ERROR("Failed to write a quad");
296 return -1;
297 }
298 qf->numquads++;
299 return 0;
300 }
301
quadfile_write_all_quads_to(quadfile_t * qf,FILE * fid)302 int quadfile_write_all_quads_to(quadfile_t* qf, FILE* fid) {
303 fitsbin_chunk_t* chunk = quads_chunk(qf);
304 if (fitsbin_write_items_to(chunk, qf->quadarray, quadfile_nquads(qf), fid)) {
305 ERROR("Failed to write %i quads", quadfile_nquads(qf));
306 return -1;
307 }
308 return 0;
309 }
310
quadfile_fix_header(quadfile_t * qf)311 int quadfile_fix_header(quadfile_t* qf) {
312 qfits_header* hdr;
313 fitsbin_t* fb = qf->fb;
314 fitsbin_chunk_t* chunk = quads_chunk(qf);
315
316 chunk->itemsize = qf->dimquads * sizeof(uint32_t);
317 chunk->nrows = qf->numquads;
318
319 hdr = fitsbin_get_primary_header(fb);
320 add_to_header(hdr, qf);
321
322 if (fitsbin_fix_primary_header(fb) ||
323 fitsbin_fix_chunk_header(fb, chunk)) {
324 ERROR("Failed to fix quad header");
325 return -1;
326 }
327 return 0;
328 }
329
quadfile_get_index_scale_upper_arcsec(const quadfile_t * qf)330 double quadfile_get_index_scale_upper_arcsec(const quadfile_t* qf) {
331 return rad2arcsec(qf->index_scale_upper);
332 }
333
quadfile_get_index_scale_lower_arcsec(const quadfile_t * qf)334 double quadfile_get_index_scale_lower_arcsec(const quadfile_t* qf) {
335 return rad2arcsec(qf->index_scale_lower);
336 }
337
quadfile_get_stars(const quadfile_t * qf,unsigned int quadid,unsigned int * stars)338 int quadfile_get_stars(const quadfile_t* qf, unsigned int quadid, unsigned int* stars) {
339 int i;
340 if (quadid >= qf->numquads) {
341 ERROR("Requested quadid %i, but number of quads is %i", quadid, qf->numquads);
342 assert(quadid < qf->numquads);
343 return -1;
344 }
345
346 for (i=0; i<qf->dimquads; i++) {
347 stars[i] = qf->quadarray[quadid * qf->dimquads + i];
348 }
349 return 0;
350 }
351
352
353
354