1 /*
2 # This file is part of libkd.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5
6 #include <unistd.h>
7 #include <stdlib.h>
8 #include <stdio.h>
9 #include <math.h>
10 #include <string.h>
11
12 #include "kdtree.h"
13 #include "kdtree_fits_io.h"
14 #include "ioutils.h"
15 #include "fitsioutils.h"
16 #include "errors.h"
17 #include "anqfits.h"
18
printHelp(char * progname)19 void printHelp(char* progname) {
20 printf("\nUsage: %s <input> <output>\n"
21 "\n", progname);
22 }
23
24
25 const char* OPTIONS = "hvq";
26
main(int argc,char ** args)27 int main(int argc, char** args) {
28 int argchar;
29 char* progname = args[0];
30 kdtree_t* kd;
31 char* infn;
32 char* outfn;
33 qfits_header* hdr;
34 qfits_header* outhdr;
35 int i, Next;
36 FILE* fout;
37 FILE* fin;
38 anbool verbose = FALSE;
39 char* err;
40 anbool force_quad = FALSE;
41 anqfits_t* anq = NULL;
42
43 while ((argchar = getopt(argc, args, OPTIONS)) != -1)
44 switch (argchar) {
45 case 'v':
46 verbose = TRUE;
47 break;
48 case 'q':
49 force_quad = TRUE;
50 break;
51 case 'h':
52 printHelp(progname);
53 exit(-1);
54 }
55
56 if (optind != argc - 2) {
57 printHelp(progname);
58 exit(-1);
59 }
60
61 infn = args[optind];
62 outfn = args[optind+1];
63
64 if (!strcmp(infn, outfn)) {
65 printf("Sorry, in-place modification of files is not supported.\n");
66 exit(-1);
67 }
68
69 if (!force_quad && ends_with(infn, ".quad.fits")) {
70 printf("\nYou don't need to fix .quad.fits files.\n"
71 " (use the -q option to try anyway.)\n");
72 exit(1);
73 }
74
75 printf("Reading kdtree from file %s ...\n", infn);
76
77 errors_start_logging_to_string();
78 kd = kdtree_fits_read(infn, NULL, &hdr);
79 err = errors_stop_logging_to_string("\n ");
80 if (!kd) {
81 printf("Failed to read kdtree from file %s:\n", infn);
82 printf(" %s\n", err);
83 free(err);
84 exit(-1);
85 }
86 free(err);
87
88 if (!kdtree_has_old_bb(kd)) {
89 printf("Kdtree %s has the correct number of bounding boxes; it doesn't need fixing.\n", infn);
90 exit(1);
91 }
92
93 if (verbose) {
94 printf("Tree name: %s\n", kd->name);
95 printf("Treetype: 0x%x\n", kd->treetype);
96 printf("Data type: %s\n", kdtree_kdtype_to_string(kdtree_datatype(kd)));
97 printf("Tree type: %s\n", kdtree_kdtype_to_string(kdtree_treetype(kd)));
98 printf("External type: %s\n", kdtree_kdtype_to_string(kdtree_exttype(kd)));
99 printf("N data points: %i\n", kd->ndata);
100 printf("Dimensions: %i\n", kd->ndim);
101 printf("Nodes: %i\n", kd->nnodes);
102 printf("Leaf nodes: %i\n", kd->nbottom);
103 printf("Non-leaf nodes: %i\n", kd->ninterior);
104 printf("Tree levels: %i\n", kd->nlevels);
105 printf("LR array: %s\n", (kd->lr ? "yes" : "no"));
106 printf("Perm array: %s\n", (kd->perm ? "yes" : "no"));
107 printf("Bounding box: %s\n", (kd->bb.any ? "yes" : "no"));
108 printf("Split plane: %s\n", (kd->split.any ? "yes" : "no"));
109 printf("Split dim: %s\n", (kd->splitdim ? "yes" : "no"));
110 printf("Data: %s\n", (kd->data.any ? "yes" : "no"));
111
112 if (kd->minval && kd->maxval) {
113 int d;
114 printf("Data ranges:\n");
115 for (d=0; d<kd->ndim; d++)
116 printf(" %i: [%g, %g]\n", d, kd->minval[d], kd->maxval[d]);
117 }
118 }
119
120 if (verbose)
121 printf("Computing bounding boxes...\n");
122 kdtree_fix_bounding_boxes(kd);
123
124 if (verbose)
125 printf("Running kdtree_check...\n");
126 if (kdtree_check(kd)) {
127 printf("kdtree_check failed.\n");
128 exit(-1);
129 }
130
131 outhdr = qfits_header_new();
132 fits_append_long_comment(outhdr, "This file was processed by the fix-bb "
133 "program, part of the Astrometry.net suite. The "
134 "extra FITS headers in the original file are "
135 "given below:");
136 fits_append_long_comment(outhdr, "---------------------------------");
137
138 for (i=0; i<qfits_header_n(hdr); i++) {
139 char key[FITS_LINESZ+1];
140 char val[FITS_LINESZ+1];
141 char com[FITS_LINESZ+1];
142 qfits_header_getitem(hdr, i, key, val, com, NULL);
143 if (!(fits_is_primary_header(key) ||
144 fits_is_table_header(key))) {
145 qfits_header_append(outhdr, key, val, com, NULL);
146 }
147 }
148 fits_append_long_comment(outhdr, "---------------------------------");
149
150 if (kdtree_fits_write(kd, outfn, outhdr)) {
151 ERROR("Failed to write output");
152 exit(-1);
153 }
154
155 if (verbose)
156 printf("Finding extra extensions...\n");
157
158 fin = fopen(infn, "rb");
159 if (!fin) {
160 SYSERROR("Failed to re-open input file %s for reading", infn);
161 exit(-1);
162 }
163 fout = fopen(outfn, "ab");
164 if (!fout) {
165 SYSERROR("Failed to re-open output file %s for writing", outfn);
166 exit(-1);
167 }
168
169 anq = anqfits_open(infn);
170 if (!anq) {
171 ERROR("Failed to open input file %s for reading", infn);
172 exit(-1);
173 }
174 Next = anqfits_n_ext(anq);
175
176 for (i=0; i<Next; i++) {
177 int hoffset, hlength;
178 int doffset, dlength;
179 int ext = i+1;
180
181 if (anqfits_is_table(anq, ext)) {
182 qfits_table* table;
183 table = anqfits_get_table(anq, ext);
184 if (table &&
185 (table->nc == 1) &&
186 kdtree_fits_column_is_kdtree(table->col[0].tlabel))
187 continue;
188 }
189 if (verbose)
190 printf("Extension %i is not part of the kdtree. Copying it verbatim.\n", ext);
191
192 hoffset = anqfits_header_start(anq, i);
193 hlength = anqfits_header_size (anq, i);
194 doffset = anqfits_data_start(anq, i);
195 dlength = anqfits_data_size (anq, i);
196
197 if (pipe_file_offset(fin, hoffset, hlength, fout) ||
198 pipe_file_offset(fin, doffset, dlength, fout)) {
199 ERROR("Failed to write extension %i verbatim", ext);
200 exit(-1);
201 }
202 }
203 fclose(fin);
204 if (fclose(fout)) {
205 SYSERROR("Failed to close output file %s", outfn);
206 exit(-1);
207 }
208 anqfits_close(anq);
209 kdtree_fits_close(kd);
210 errors_free();
211
212 printf("Fixed file %s was written successfully.\n", outfn);
213
214 return 0;
215 }
216