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 "starkd.h"
24 #include "boilerplate.h"
25 #include "errors.h"
26 #include "log.h"
27 #include "quad-utils.h"
28 #include "quad-builder.h"
29 #include "allquads.h"
30 #include "starkd.h"
31 #include "startree.h"
32 #include "kdtree.h"
33 #include "codekd.h"
34 #include "codetree.h"
35 #include "unpermute-quads.h"
36 #include "unpermute-stars.h"
37 #include "index.h"
38 #include "merge-index.h"
39 
40 //#include "build-index.h"
41 
42 const char* OPTIONS = "hi:o:u:l:d:I:v";
43 
print_help(char * progname)44 static void print_help(char* progname) {
45     BOILERPLATE_HELP_HEADER(stdout);
46     printf("\nUsage: %s\n"
47            "      -i <input-filename>    (FITS reference catalog)\n"
48            "      -o <index-filename>    (Index file output filename)\n"
49            "     [-u <scale>]    upper bound of quad scale (arcmin)\n"
50            "     [-l <scale>]    lower bound of quad scale (arcmin)\n"
51            "     [-d <dimquads>] number of stars in a \"quad\".\n"
52            "     [-I <unique-id>] set the unique ID of this index\n\n"
53            "\nReads skdt, writes {code, quad}.\n\n"
54            , progname);
55 }
56 
57 
main(int argc,char ** argv)58 int main(int argc, char** argv) {
59     int argchar;
60     allquads_t* aq;
61     int loglvl = LOG_MSG;
62     int i, N;
63     char* catfn = NULL;
64 
65     startree_t* starkd;
66     fitstable_t* cat;
67 
68     char* racol = NULL;
69     char* deccol = NULL;
70     int datatype = KDT_DATA_DOUBLE;
71     int treetype = KDT_TREE_DOUBLE;
72     int buildopts = 0;
73     int Nleaf = 0;
74 
75     char* indexfn = NULL;
76     //index_params_t* p;
77 
78     aq = allquads_init();
79     aq->skdtfn = "allquads.skdt";
80     aq->codefn = "allquads.code";
81     aq->quadfn = "allquads.quad";
82 
83     while ((argchar = getopt (argc, argv, OPTIONS)) != -1)
84         switch (argchar) {
85         case 'v':
86             loglvl++;
87             break;
88         case 'd':
89             aq->dimquads = atoi(optarg);
90             break;
91         case 'I':
92             aq->id = atoi(optarg);
93             break;
94         case 'h':
95             print_help(argv[0]);
96             exit(0);
97         case 'i':
98             catfn = optarg;
99             break;
100         case 'o':
101             indexfn = optarg;
102             break;
103         case 'u':
104             aq->quad_d2_upper = arcmin2distsq(atof(optarg));
105             aq->use_d2_upper = TRUE;
106             break;
107         case 'l':
108             aq->quad_d2_lower = arcmin2distsq(atof(optarg));
109             aq->use_d2_lower = TRUE;
110             break;
111         default:
112             return -1;
113         }
114 
115     log_init(loglvl);
116 
117     if (!catfn || !indexfn) {
118         printf("Specify in & out filenames, bonehead!\n");
119         print_help(argv[0]);
120         exit( -1);
121     }
122 
123     if (optind != argc) {
124         print_help(argv[0]);
125         printf("\nExtra command-line args were given: ");
126         for (i=optind; i<argc; i++) {
127             printf("%s ", argv[i]);
128         }
129         printf("\n");
130         exit(-1);
131     }
132 
133     if (!aq->id)
134         logmsg("Warning: you should set the unique-id for this index (with -I).\n");
135 
136     if (aq->dimquads > DQMAX) {
137         ERROR("Quad dimension %i exceeds compiled-in max %i.\n", aq->dimquads, DQMAX);
138         exit(-1);
139     }
140     aq->dimcodes = dimquad2dimcode(aq->dimquads);
141 
142     // Read reference catalog, write star kd-tree
143     logmsg("Building star kdtree: reading %s, writing to %s\n", catfn, aq->skdtfn);
144 
145     logverb("Reading star catalogue...");
146     cat = fitstable_open(catfn);
147     if (!cat) {
148         ERROR("Couldn't read catalog");
149         exit(-1);
150     }
151     logmsg("Got %i stars\n", fitstable_nrows(cat));
152     starkd = startree_build(cat, racol, deccol, datatype, treetype,
153                             buildopts, Nleaf, argv, argc);
154     if (!starkd) {
155         ERROR("Failed to create star kdtree");
156         exit(-1);
157     }
158 
159     logmsg("Star kd-tree contains %i data points in dimension %i\n",
160            startree_N(starkd), startree_D(starkd));
161     N = startree_N(starkd);
162     for (i=0; i<N; i++) {
163         double ra,dec;
164         int ok;
165         ok = startree_get_radec(starkd, i, &ra, &dec);
166         logmsg("  data %i: ok %i, RA,Dec %g, %g\n", i, ok, ra, dec);
167     }
168 
169     if (startree_write_to_file(starkd, aq->skdtfn)) {
170         ERROR("Failed to write star kdtree");
171         exit(-1);
172     }
173     startree_close(starkd);
174     fitstable_close(cat);
175     logmsg("Wrote star kdtree to %s\n", aq->skdtfn);
176 
177     logmsg("Running allquads...\n");
178     if (allquads_open_outputs(aq)) {
179         exit(-1);
180     }
181 
182     if (allquads_create_quads(aq)) {
183         exit(-1);
184     }
185 
186     if (allquads_close(aq)) {
187         exit(-1);
188     }
189 
190     logmsg("allquads: wrote %s, %s\n", aq->quadfn, aq->codefn);
191 
192     // build-index:
193     //build_index_defaults(&p);
194 
195     // codetree
196     /*
197      if (step_codetree(p, codes, &codekd,
198      codefn, &ckdtfn, tempfiles))
199      return -1;
200      */
201     char* ckdtfn=NULL;
202     char* tempdir = NULL;
203 
204     ckdtfn = create_temp_file("ckdt", tempdir);
205     logmsg("Creating code kdtree, reading %s, writing to %s\n", aq->codefn, ckdtfn);
206     if (codetree_files(aq->codefn, ckdtfn, 0, 0, 0, 0, argv, argc)) {
207         ERROR("codetree failed");
208         return -1;
209     }
210 
211     char* skdt2fn=NULL;
212     char* quad2fn=NULL;
213 
214     // unpermute-stars
215     logmsg("Unpermute-stars...\n");
216     skdt2fn = create_temp_file("skdt2", tempdir);
217     quad2fn = create_temp_file("quad2", tempdir);
218 
219     logmsg("Unpermuting stars from %s and %s to %s and %s\n",
220            aq->skdtfn, aq->quadfn, skdt2fn, quad2fn);
221     if (unpermute_stars_files(aq->skdtfn, aq->quadfn, skdt2fn, quad2fn,
222                               TRUE, FALSE, argv, argc)) {
223         ERROR("Failed to unpermute-stars");
224         return -1;
225     }
226 
227     allquads_free(aq);
228 
229     // unpermute-quads
230     /*
231      if (step_unpermute_quads(p, quads2, codekd, &quads3, &codekd2,
232      quad2fn, ckdtfn, &quad3fn, &ckdt2fn, tempfiles))
233      return -1;
234      */
235     char* quad3fn=NULL;
236     char* ckdt2fn=NULL;
237 
238     ckdt2fn = create_temp_file("ckdt2", tempdir);
239     quad3fn = create_temp_file("quad3", tempdir);
240     logmsg("Unpermuting quads from %s and %s to %s and %s\n",
241            quad2fn, ckdtfn, quad3fn, ckdt2fn);
242     if (unpermute_quads_files(quad2fn, ckdtfn,
243                               quad3fn, ckdt2fn, argv, argc)) {
244         ERROR("Failed to unpermute-quads");
245         return -1;
246     }
247 
248     // index
249     /*
250      if (step_merge_index(p, codekd2, quads3, starkd2, p_index,
251      ckdt2fn, quad3fn, skdt2fn, indexfn))
252      return -1;
253      */
254     quadfile_t* quad;
255     codetree_t* code;
256     startree_t* star;
257 
258     logmsg("Merging %s and %s and %s to %s\n",
259            quad3fn, ckdt2fn, skdt2fn, indexfn);
260     if (merge_index_open_files(quad3fn, ckdt2fn, skdt2fn,
261                                &quad, &code, &star)) {
262         ERROR("Failed to open index files for merging");
263         return -1;
264     }
265     if (merge_index(quad, code, star, indexfn)) {
266         ERROR("Failed to write merged index");
267         return -1;
268     }
269     codetree_close(code);
270     startree_close(star);
271     quadfile_close(quad);
272 
273     printf("Done.\n");
274 
275     free(ckdtfn);
276     free(skdt2fn);
277     free(quad2fn);
278     free(ckdt2fn);
279     free(quad3fn);
280 
281     return 0;
282 }
283 
284