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