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 <math.h>
7 #include <stdio.h>
8 #include <errno.h>
9 #include <string.h>
10 #include <unistd.h>
11 #include <stdint.h>
12 #include <limits.h>
13 #include <sys/time.h>
14 #include <sys/resource.h>
15 #include <assert.h>
16
17 #include "starutil.h"
18 #include "codefile.h"
19 #include "mathutil.h"
20 #include "quadfile.h"
21 #include "kdtree.h"
22 #include "fitsioutils.h"
23 #include "anqfits.h"
24 #include "starkd.h"
25 #include "boilerplate.h"
26 #include "errors.h"
27 #include "log.h"
28 #include "quad-utils.h"
29 #include "quad-builder.h"
30 #include "allquads.h"
31
allquads_init()32 allquads_t* allquads_init() {
33 allquads_t* aq = calloc(1, sizeof(allquads_t));
34 aq->dimquads = 4;
35 return aq;
36 }
37
allquads_free(allquads_t * aq)38 void allquads_free(allquads_t* aq) {
39 free(aq);
40 }
41
add_quad(quadbuilder_t * qb,unsigned int * quad,void * token)42 static void add_quad(quadbuilder_t* qb, unsigned int* quad, void* token) {
43 allquads_t* aq = token;
44 if (log_get_level() > LOG_VERB) {
45 int k;
46 debug("quad: ");
47 for (k=0; k<qb->dimquads; k++)
48 debug("%-6i ", quad[k]);
49 logverb("\n");
50 }
51 quad_write_const(aq->codes, aq->quads, quad, aq->starkd, qb->dimquads, aq->dimcodes);
52 }
53
check_AB(quadbuilder_t * qb,pquad_t * pq,void * token)54 static anbool check_AB(quadbuilder_t* qb, pquad_t* pq, void* token) {
55 allquads_t* aq = token;
56 debug("check_AB: iA=%i, iB=%i, aq->starA=%i\n", pq->iA, pq->iB, aq->starA);
57 return (pq->iA == aq->starA);
58 }
59
60 /*
61 #include "codefile.h"
62 static bool check_full_quad(quadbuilder_t* qb, unsigned int* quad, int nstars, void* token) {
63 allquads_t* aq = token;
64 double code[4];
65 double xyz[12];
66 double R1, R2;
67 startree_get(aq->starkd, quad[0], xyz+0);
68 startree_get(aq->starkd, quad[1], xyz+3);
69 startree_get(aq->starkd, quad[2], xyz+6);
70 startree_get(aq->starkd, quad[3], xyz+9);
71 codefile_compute_star_code(xyz, code, 4);
72 logmsg("Code: %g, %g, %g, %g\n", code[0], code[1], code[2], code[3]);
73 R1 = hypot(code[0] - 0.5, code[1] - 0.5);
74 R2 = hypot(code[2] - 0.5, code[3] - 0.5);
75 logmsg("Radii %g, %g\n", R1, R2);
76 return TRUE;
77 }
78 */
79
allquads_create_quads(allquads_t * aq)80 int allquads_create_quads(allquads_t* aq) {
81 quadbuilder_t* qb;
82 int i, N;
83 double* xyz;
84
85 qb = quadbuilder_init();
86
87 qb->quadd2_high = aq->quad_d2_upper;
88 qb->quadd2_low = aq->quad_d2_lower;
89 qb->check_scale_high = aq->use_d2_upper;
90 qb->check_scale_low = aq->use_d2_lower;
91
92 qb->dimquads = aq->dimquads;
93 qb->add_quad = add_quad;
94 qb->add_quad_token = aq;
95
96 N = startree_N(aq->starkd);
97
98 logverb("check scale high: %i\n", qb->check_scale_high);
99
100 if (!qb->check_scale_high) {
101 int* inds;
102 inds = malloc(N * sizeof(int));
103 for (i=0; i<N; i++)
104 inds[i] = i;
105 xyz = malloc(3 * N * sizeof(double));
106 kdtree_copy_data_double(aq->starkd->tree, 0, N, xyz);
107
108 qb->starxyz = xyz;
109 qb->starinds = inds;
110 qb->Nstars = N;
111
112 quadbuilder_create(qb);
113
114 free(xyz);
115 free(inds);
116 } else {
117 int nq;
118 int lastgrass = 0;
119
120 /*
121 xyz = malloc(3 * N * sizeof(double));
122 kdtree_copy_data_double(aq->starkd->tree, 0, N, xyz);
123 */
124
125 // star A = i
126 nq = aq->quads->numquads;
127 for (i=0; i<N; i++) {
128 double xyzA[3];
129 int* inds;
130 int NR;
131
132 int grass = (i*80 / N);
133 if (grass != lastgrass) {
134 printf(".");
135 fflush(stdout);
136 lastgrass = grass;
137 }
138
139 startree_get(aq->starkd, i, xyzA);
140
141 startree_search_for(aq->starkd, xyzA, aq->quad_d2_upper,
142 &xyz, NULL, &inds, &NR);
143
144 /*
145 startree_search_for(aq->starkd, xyzA, aq->quad_d2_upper,
146 NULL, NULL, &inds, &NR);
147 */
148
149 logverb("Star %i of %i: found %i stars in range\n", i+1, N, NR);
150 aq->starA = i;
151 qb->starxyz = xyz;
152 qb->starinds = inds;
153 qb->Nstars = NR;
154 qb->check_AB_stars = check_AB;
155 qb->check_AB_stars_token = aq;
156 //qb->check_full_quad = check_full_quad;
157 //qb->check_full_quad_token = aq;
158
159 quadbuilder_create(qb);
160
161 logverb("Star %i of %i: wrote %i quads for this star, total %i so far.\n", i+1, N, aq->quads->numquads - nq, aq->quads->numquads);
162 free(inds);
163 free(xyz);
164 }
165 //
166 //free(xyz);
167
168 printf("\n");
169 }
170
171 quadbuilder_free(qb);
172 return 0;
173 }
174
allquads_close(allquads_t * aq)175 int allquads_close(allquads_t* aq) {
176 startree_close(aq->starkd);
177
178 // fix output file headers.
179 if (quadfile_fix_header(aq->quads) ||
180 quadfile_close(aq->quads)) {
181 ERROR("Couldn't write quad output file");
182 return -1;
183 }
184 if (codefile_fix_header(aq->codes) ||
185 codefile_close(aq->codes)) {
186 ERROR("Couldn't write code output file");
187 return -1;
188 }
189 return 0;
190 }
191
allquads_open_outputs(allquads_t * aq)192 int allquads_open_outputs(allquads_t* aq) {
193 int hp, hpnside;
194 qfits_header* hdr;
195
196 printf("Reading star kdtree %s ...\n", aq->skdtfn);
197 aq->starkd = startree_open(aq->skdtfn);
198 if (!aq->starkd) {
199 ERROR("Failed to open star kdtree %s\n", aq->skdtfn);
200 return -1;
201 }
202 printf("Star tree contains %i objects.\n", startree_N(aq->starkd));
203
204 printf("Will write to quad file %s and code file %s\n", aq->quadfn, aq->codefn);
205 aq->quads = quadfile_open_for_writing(aq->quadfn);
206 if (!aq->quads) {
207 ERROR("Couldn't open file %s to write quads.\n", aq->quadfn);
208 return -1;
209 }
210 aq->codes = codefile_open_for_writing(aq->codefn);
211 if (!aq->codes) {
212 ERROR("Couldn't open file %s to write codes.\n", aq->quadfn);
213 return -1;
214 }
215
216 aq->quads->dimquads = aq->dimquads;
217 aq->codes->dimcodes = aq->dimcodes;
218
219 if (aq->id) {
220 aq->quads->indexid = aq->id;
221 aq->codes->indexid = aq->id;
222 }
223
224 // get the "HEALPIX" header from the skdt and put it in the code and quad headers.
225 hp = qfits_header_getint(startree_header(aq->starkd), "HEALPIX", -1);
226 if (hp == -1) {
227 logmsg("Warning: skdt does not contain \"HEALPIX\" header. Code and quad files will not contain this header either.\n");
228 }
229 aq->quads->healpix = hp;
230 aq->codes->healpix = hp;
231 // likewise "HPNSIDE"
232 hpnside = qfits_header_getint(startree_header(aq->starkd), "HPNSIDE", 1);
233 aq->quads->hpnside = hpnside;
234 aq->codes->hpnside = hpnside;
235
236 hdr = quadfile_get_header(aq->quads);
237 qfits_header_add(hdr, "CXDX", "T", "All codes have the property cx<=dx.", NULL);
238 qfits_header_add(hdr, "CXDXLT1", "T", "All codes have the property cx+dx<=1.", NULL);
239 qfits_header_add(hdr, "MIDHALF", "T", "All codes have the property cx+dx<=1.", NULL);
240 qfits_header_add(hdr, "CIRCLE", "T", "Codes live in the circle, not the box.", NULL);
241
242 hdr = codefile_get_header(aq->codes);
243 qfits_header_add(hdr, "CXDX", "T", "All codes have the property cx<=dx.", NULL);
244 qfits_header_add(hdr, "CXDXLT1", "T", "All codes have the property cx+dx<=1.", NULL);
245 qfits_header_add(hdr, "MIDHALF", "T", "All codes have the property cx+dx<=1.", NULL);
246 qfits_header_add(hdr, "CIRCLE", "T", "Codes live in the circle, not the box.", NULL);
247
248 if (quadfile_write_header(aq->quads)) {
249 ERROR("Couldn't write headers to quads file %s\n", aq->quadfn);
250 return -1;
251 }
252 if (codefile_write_header(aq->codes)) {
253 ERROR("Couldn't write headers to code file %s\n", aq->codefn);
254 return -1;
255 }
256
257 if (!aq->use_d2_lower)
258 aq->quad_d2_lower = 0.0;
259 if (!aq->use_d2_upper)
260 aq->quad_d2_upper = 10.0;
261
262 aq->codes->numstars = startree_N(aq->starkd);
263 aq->codes->index_scale_upper = distsq2rad(aq->quad_d2_upper);
264 aq->codes->index_scale_lower = distsq2rad(aq->quad_d2_lower);
265
266 aq->quads->numstars = aq->codes->numstars;
267 aq->quads->index_scale_upper = aq->codes->index_scale_upper;
268 aq->quads->index_scale_lower = aq->codes->index_scale_lower;
269 return 0;
270 }
271
272