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