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 <string.h>
6 #include <math.h>
7 #include <assert.h>
8
9 #include "plotindex.h"
10 #include "cairoutils.h"
11 #include "ioutils.h"
12 #include "log.h"
13 #include "errors.h"
14 #include "starutil.h"
15 #include "index.h"
16 #include "qidxfile.h"
17 #include "permutedsort.h"
18
19 DEFINE_PLOTTER(index);
20
plot_index_get(plot_args_t * pargs)21 plotindex_t* plot_index_get(plot_args_t* pargs) {
22 return plotstuff_get_config(pargs, "index");
23 }
24
plot_index_init(plot_args_t * plotargs)25 void* plot_index_init(plot_args_t* plotargs) {
26 plotindex_t* args = calloc(1, sizeof(plotindex_t));
27 args->indexes = pl_new(16);
28 args->qidxes = pl_new(16);
29 args->stars = TRUE;
30 args->quads = TRUE;
31 args->fill = FALSE;
32 return args;
33 }
34
pad_qidxes(plotindex_t * args)35 static void pad_qidxes(plotindex_t* args) {
36 while ((pl_size(args->qidxes)) < pl_size(args->indexes))
37 pl_append(args->qidxes, NULL);
38 }
39
plot_quad_xy(cairo_t * cairo,double * quadxy,int dimquads)40 void plot_quad_xy(cairo_t* cairo, double* quadxy, int dimquads) {
41 int k;
42 double cx, cy;
43 double theta[DQMAX];
44 int* perm;
45
46 cx = cy = 0.0;
47 for (k=0; k<dimquads; k++) {
48 cx += quadxy[2*k+0];
49 cy += quadxy[2*k+1];
50 }
51 cx /= dimquads;
52 cy /= dimquads;
53
54 for (k=0; k<dimquads; k++)
55 theta[k] = atan2(quadxy[2*k+1] - cy, quadxy[2*k+0] - cx);
56 perm = permuted_sort(theta, sizeof(double), compare_doubles_asc, NULL, dimquads);
57 for (k=0; k<dimquads; k++) {
58 double px,py;
59 px = quadxy[2 * perm[k] + 0];
60 py = quadxy[2 * perm[k] + 1];
61 if (k == 0) {
62 cairo_move_to(cairo, px, py);
63 } else {
64 cairo_line_to(cairo, px, py);
65 }
66 }
67 free(perm);
68 cairo_close_path(cairo);
69 }
70
plotquad(cairo_t * cairo,plot_args_t * pargs,plotindex_t * args,index_t * index,int quadnum,int DQ)71 static void plotquad(cairo_t* cairo, plot_args_t* pargs, plotindex_t* args, index_t* index, int quadnum, int DQ) {
72 int k;
73 unsigned int stars[DQMAX];
74 double ra, dec;
75 double px, py;
76 double xy[DQMAX*2];
77 int N;
78
79 quadfile_get_stars(index->quads, quadnum, stars);
80 N = 0;
81 for (k=0; k<DQ; k++) {
82 if (startree_get_radec(index->starkd, stars[k], &ra, &dec)) {
83 ERROR("Failed to get RA,Dec for star %i\n", stars[k]);
84 continue;
85 }
86 if (!plotstuff_radec2xy(pargs, ra, dec, &px, &py)) {
87 ERROR("Failed to convert RA,Dec %g,%g to pixels for quad %i\n", ra, dec, quadnum);
88 continue;
89 }
90 xy[2*k + 0] = px;
91 xy[2*k + 1] = py;
92 N++;
93 }
94 if (N < 3)
95 return;
96 plot_quad_xy(cairo, xy, N);
97 if (args->fill)
98 cairo_fill(cairo);
99 else
100 cairo_stroke(cairo);
101 }
102
plot_index_plotquad(cairo_t * cairo,plot_args_t * pargs,plotindex_t * args,index_t * index,int quadnum,int DQ)103 void plot_index_plotquad(cairo_t* cairo, plot_args_t* pargs, plotindex_t* args, index_t* index, int quadnum, int DQ) {
104 plotquad(cairo, pargs, args, index, quadnum, DQ);
105 }
106
plot_index_plot(const char * command,cairo_t * cairo,plot_args_t * pargs,void * baton)107 int plot_index_plot(const char* command,
108 cairo_t* cairo, plot_args_t* pargs, void* baton) {
109 plotindex_t* args = (plotindex_t*)baton;
110 int i;
111 double ra, dec, radius;
112 double xyz[3];
113 double r2;
114
115 pad_qidxes(args);
116
117 plotstuff_builtin_apply(cairo, pargs);
118
119 if (plotstuff_get_radec_center_and_radius(pargs, &ra, &dec, &radius)) {
120 ERROR("Failed to get RA,Dec center and radius");
121 return -1;
122 }
123 radecdeg2xyzarr(ra, dec, xyz);
124 r2 = deg2distsq(radius);
125 logmsg("Field RA,Dec,radius = (%g,%g), %g deg\n", ra, dec, radius);
126 logmsg("distsq: %g\n", r2);
127
128 for (i=0; i<pl_size(args->indexes); i++) {
129 index_t* index = pl_get(args->indexes, i);
130 int j, N;
131 int DQ;
132 double px,py;
133
134 if (args->stars) {
135 // plot stars
136 double* radecs = NULL;
137 startree_search_for(index->starkd, xyz, r2, NULL, &radecs, NULL, &N);
138 if (N) {
139 assert(radecs);
140 }
141 logmsg("Found %i stars in range in index %s\n", N, index->indexname);
142 for (j=0; j<N; j++) {
143 logverb(" RA,Dec (%g,%g) -> x,y (%g,%g)\n", radecs[2*j], radecs[2*j+1], px, py);
144 if (!plotstuff_radec2xy(pargs, radecs[j*2], radecs[j*2+1], &px, &py)) {
145 ERROR("Failed to convert RA,Dec %g,%g to pixels\n", radecs[j*2], radecs[j*2+1]);
146 continue;
147 }
148 cairoutils_draw_marker(cairo, pargs->marker, px, py, pargs->markersize);
149 cairo_stroke(cairo);
150 }
151 free(radecs);
152 }
153 if (args->quads) {
154 DQ = index_get_quad_dim(index);
155 qidxfile* qidx = pl_get(args->qidxes, i);
156 if (qidx) {
157 int* stars;
158 int Nstars;
159 il* quadlist = il_new(256);
160
161 // find stars in range.
162 startree_search_for(index->starkd, xyz, r2, NULL, NULL, &stars, &Nstars);
163 logmsg("Found %i stars in range of index %s\n", N, index->indexname);
164 logmsg("Using qidx file.\n");
165 // find quads that each star is a member of.
166 for (j=0; j<Nstars; j++) {
167 uint32_t* quads;
168 int Nquads;
169 int k;
170 if (qidxfile_get_quads(qidx, stars[j], &quads, &Nquads)) {
171 ERROR("Failed to get quads for star %i\n", stars[j]);
172 return -1;
173 }
174 for (k=0; k<Nquads; k++)
175 il_insert_unique_ascending(quadlist, quads[k]);
176 }
177 for (j=0; j<il_size(quadlist); j++) {
178 plotquad(cairo, pargs, args, index, il_get(quadlist, j), DQ);
179 }
180
181 } else {
182 // plot quads
183 N = index_nquads(index);
184 for (j=0; j<N; j++) {
185 plotquad(cairo, pargs, args, index, j, DQ);
186 }
187 }
188 }
189 }
190 return 0;
191 }
192
plot_index_add_qidx_file(plotindex_t * args,const char * fn)193 int plot_index_add_qidx_file(plotindex_t* args, const char* fn) {
194 int i;
195 qidxfile* qidx = qidxfile_open(fn);
196 if (!qidx) {
197 ERROR("Failed to open quad index file \"%s\"", fn);
198 return -1;
199 }
200 pad_qidxes(args);
201 i = pl_size(args->indexes) - 1;
202 pl_set(args->qidxes, i, qidx);
203 return 0;
204 }
205
plot_index_add_file(plotindex_t * args,const char * fn)206 int plot_index_add_file(plotindex_t* args, const char* fn) {
207 index_t* index = index_load(fn, 0, NULL);
208 if (!index) {
209 ERROR("Failed to open index \"%s\"", fn);
210 return -1;
211 }
212 pl_append(args->indexes, index);
213 return 0;
214 }
215
plot_index_command(const char * cmd,const char * cmdargs,plot_args_t * pargs,void * baton)216 int plot_index_command(const char* cmd, const char* cmdargs,
217 plot_args_t* pargs, void* baton) {
218 plotindex_t* args = (plotindex_t*)baton;
219 if (streq(cmd, "index_file")) {
220 const char* fn = cmdargs;
221 return plot_index_add_file(args, fn);
222 } else if (streq(cmd, "index_qidxfile")) {
223 const char* fn = cmdargs;
224 return plot_index_add_qidx_file(args, fn);
225 } else if (streq(cmd, "index_draw_stars")) {
226 args->stars = atoi(cmdargs);
227 } else if (streq(cmd, "index_draw_quads")) {
228 args->quads = atoi(cmdargs);
229 } else if (streq(cmd, "index_fill")) {
230 args->fill = atoi(cmdargs);
231 } else {
232 ERROR("Did not understand command \"%s\"", cmd);
233 return -1;
234 }
235 return 0;
236 }
237
plot_index_free(plot_args_t * plotargs,void * baton)238 void plot_index_free(plot_args_t* plotargs, void* baton) {
239 plotindex_t* args = (plotindex_t*)baton;
240 int i;
241 for (i=0; i<pl_size(args->indexes); i++) {
242 index_t* index = pl_get(args->indexes, i);
243 index_free(index);
244 }
245 pl_free(args->indexes);
246 for (i=0; i<pl_size(args->qidxes); i++) {
247 qidxfile* qidx = pl_get(args->qidxes, i);
248 qidxfile_close(qidx);
249 }
250 pl_free(args->qidxes);
251 free(args);
252 }
253
254