1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5 #include <assert.h>
6 #include "startree.h"
7 #include "kdtree.h"
8 #include "errors.h"
9 #include "log.h"
10 #include "starutil.h"
11 #include "an-bool.h"
12 #include "fitsioutils.h"
13 #include "boilerplate.h"
14 #include "fitstable.h"
15
startree_has_tagalong_data(const fitstable_t * intab)16 anbool startree_has_tagalong_data(const fitstable_t* intab) {
17 // don't include RA,Dec.
18 return fitstable_n_fits_columns(intab) > 2;
19 }
20
startree_write_tagalong_table(fitstable_t * intab,fitstable_t * outtab,const char * racol,const char * deccol,int * indices,anbool remove_radec_columns)21 int startree_write_tagalong_table(fitstable_t* intab, fitstable_t* outtab,
22 const char* racol, const char* deccol,
23 int* indices,
24 anbool remove_radec_columns) {
25 int i, R, NB, N;
26 qfits_header* hdr;
27
28 fitstable_clear_table(intab);
29 fitstable_add_fits_columns_as_struct(intab);
30 fitstable_copy_columns(intab, outtab);
31
32 if (remove_radec_columns) {
33 if (!racol)
34 racol = "RA";
35 if (!deccol)
36 deccol = "DEC";
37 fitstable_remove_column(outtab, racol);
38 fitstable_remove_column(outtab, deccol);
39 }
40 fitstable_read_extension(intab, 1);
41 hdr = fitstable_get_header(outtab);
42 qfits_header_add(hdr, "AN_FILE", AN_FILETYPE_TAGALONG, "Extra data for stars", NULL);
43 if (fitstable_write_header(outtab)) {
44 ERROR("Failed to write tag-along data header");
45 return -1;
46 }
47 N = fitstable_nrows(intab);
48 R = fitstable_row_size(intab);
49
50 if (indices) {
51 if (!remove_radec_columns) {
52 // row-by-row raw data copy; read whole tag-along array into memory
53 char* data = malloc((size_t)N * (size_t)R);
54 // FIXME -- could read row-by-row if the malloc fails.....
55 if (!data) {
56 ERROR("Failed to allocate enough memory to read full tag-along table");
57 return -1;
58 }
59 printf("Reading tag-along table...\n");
60 if (fitstable_read_nrows_data(intab, 0, N, data)) {
61 ERROR("Failed to read tag-along table");
62 free(data);
63 return -1;
64 }
65 printf("Writing tag-along table...\n");
66 for (i=0; i<N; i++) {
67 if (fitstable_write_row_data(outtab, data + (size_t)indices[i]*(size_t)R)) {
68 ERROR("Failed to write a row of data");
69 free(data);
70 return -1;
71 }
72 }
73 free(data);
74
75 } else {
76 if (fitstable_copy_rows_data(intab, indices, N, outtab)) {
77 ERROR("Failed to copy tag-along table rows from input to output");
78 return -1;
79 }
80 }
81 } else {
82 char* buf;
83
84 NB = 1000;
85 logverb("Input row size: %i, output row size: %i\n", R, fitstable_row_size(outtab));
86 buf = malloc(NB * R);
87
88 for (i=0; i<N; i+=NB) {
89 int nr = NB;
90 if (i+NB > N)
91 nr = N - i;
92 if (fitstable_read_structs(intab, buf, R, i, nr)) {
93 ERROR("Failed to read tag-along data from catalog");
94 return -1;
95 }
96 if (fitstable_write_structs(outtab, buf, R, nr)) {
97 ERROR("Failed to write tag-along data");
98 return -1;
99 }
100 }
101 free(buf);
102 }
103 if (fitstable_fix_header(outtab)) {
104 ERROR("Failed to fix tag-along data header");
105 return -1;
106 }
107 return 0;
108 }
109
startree_build(fitstable_t * intable,const char * racol,const char * deccol,int datatype,int treetype,int buildopts,int Nleaf,char ** args,int argc)110 startree_t* startree_build(fitstable_t* intable,
111 const char* racol, const char* deccol,
112 // keep RA,Dec in the tag-along table?
113 //anbool keep_radec,
114 // KDT_DATA_*, KDT_TREE_*
115 int datatype, int treetype,
116 // KD_BUILD_*
117 int buildopts,
118 int Nleaf,
119 char** args, int argc) {
120 double* ra = NULL;
121 double* dec = NULL;
122 double* xyz = NULL;
123 int N;
124 startree_t* starkd = NULL;
125 int tt;
126 int d;
127 double low[3];
128 double high[3];
129 qfits_header* hdr;
130 qfits_header* inhdr;
131 int i;
132
133 if (!racol)
134 racol = "RA";
135 if (!deccol)
136 deccol = "DEC";
137 if (!datatype)
138 datatype = KDT_DATA_U32;
139 if (!treetype)
140 treetype = KDT_TREE_U32;
141 if (!buildopts)
142 buildopts = KD_BUILD_SPLIT;
143 if (!Nleaf)
144 Nleaf = 25;
145
146 ra = fitstable_read_column(intable, racol, TFITS_BIN_TYPE_D);
147 if (!ra) {
148 ERROR("Failed to read RA from column %s", racol);
149 goto bailout;
150 }
151 dec = fitstable_read_column(intable, deccol, TFITS_BIN_TYPE_D);
152 if (!dec) {
153 ERROR("Failed to read RA from column %s", racol);
154 goto bailout;
155 }
156 printf("First RA,Dec: %g,%g\n", ra[0], dec[0]);
157 N = fitstable_nrows(intable);
158 xyz = malloc(N * 3 * sizeof(double));
159 if (!xyz) {
160 SYSERROR("Failed to malloc xyz array to build startree");
161 goto bailout;
162 }
163 radecdeg2xyzarrmany(ra, dec, xyz, N);
164 free(ra);
165 ra = NULL;
166 free(dec);
167 dec = NULL;
168 printf("First x,y,z: %g,%g,%g\n", xyz[0], xyz[1], xyz[2]);
169
170 starkd = startree_new();
171 if (!starkd) {
172 ERROR("Failed to allocate startree");
173 free(xyz);
174 goto bailout;
175 }
176 tt = kdtree_kdtypes_to_treetype(KDT_EXT_DOUBLE, treetype, datatype);
177 printf("Treetype: 0x%x\n", tt);
178 starkd->tree = kdtree_new(N, 3, Nleaf);
179 for (d=0; d<3; d++) {
180 low[d] = -1.0;
181 high[d] = 1.0;
182 }
183 kdtree_set_limits(starkd->tree, low, high);
184 logverb("Building star kdtree...\n");
185 starkd->tree = kdtree_build(starkd->tree, xyz, N, 3, Nleaf, tt, buildopts);
186 if (!starkd->tree) {
187 ERROR("Failed to build star kdtree");
188 startree_close(starkd);
189 starkd = NULL;
190 free(xyz);
191 goto bailout;
192 }
193 starkd->tree->name = strdup(STARTREE_NAME);
194
195 printf("After kdtree_build:\n");
196 kdtree_print(starkd->tree);
197 {
198 double* treed = kdtree_get_data(starkd->tree, 0);
199 printf("First data elements in tree: %g,%g,%g\n", treed[0], treed[1], treed[2]);
200 }
201
202 inhdr = fitstable_get_primary_header(intable);
203 hdr = startree_header(starkd);
204 an_fits_copy_header(inhdr, hdr, "HEALPIX");
205 an_fits_copy_header(inhdr, hdr, "HPNSIDE");
206 an_fits_copy_header(inhdr, hdr, "ALLSKY");
207 an_fits_copy_header(inhdr, hdr, "JITTER");
208 an_fits_copy_header(inhdr, hdr, "CUTNSIDE");
209 an_fits_copy_header(inhdr, hdr, "CUTMARG");
210 an_fits_copy_header(inhdr, hdr, "CUTDEDUP");
211 an_fits_copy_header(inhdr, hdr, "CUTNSWEP");
212 //fits_copy_header(inhdr, hdr, "CUTBAND");
213 //fits_copy_header(inhdr, hdr, "CUTMINMG");
214 //fits_copy_header(inhdr, hdr, "CUTMAXMG");
215 BOILERPLATE_ADD_FITS_HEADERS(hdr);
216 qfits_header_add(hdr, "HISTORY", "This file was created by the command-line:", NULL, NULL);
217 fits_add_args(hdr, args, argc);
218 qfits_header_add(hdr, "HISTORY", "(end of command line)", NULL, NULL);
219 qfits_header_add(hdr, "HISTORY", "** History entries copied from the input file:", NULL, NULL);
220 fits_copy_all_headers(inhdr, hdr, "HISTORY");
221 qfits_header_add(hdr, "HISTORY", "** End of history entries.", NULL, NULL);
222 for (i=1;; i++) {
223 char key[16];
224 int n;
225 sprintf(key, "SWEEP%i", i);
226 n = qfits_header_getint(inhdr, key, -1);
227 if (n == -1)
228 break;
229 an_fits_copy_header(inhdr, hdr, key);
230 }
231
232 bailout:
233 if (ra)
234 free(ra);
235 if (dec)
236 free(dec);
237 // NOOO don't free xyz -- it belongs to the kdtree!
238 //if (xyz)
239 //free(xyz);
240 return starkd;
241 }
242
243