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