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 <unistd.h>
6 #include <stdio.h>
7 #include <string.h>
8 #include <math.h>
9 #include <assert.h>
10
11 /**
12
13 wget "http://antwrp.gsfc.nasa.gov/apod/image/0403/cmsky_cortner_full.jpg"
14 solve-field --backend-config backend.cfg -v --keep-xylist %s.xy --continue --scale-low 10 --scale-units degwidth --no-tweak cmsky_cortner_full.jpg
15 # index-219 ==> 159 index objs.
16 cp cmsky_cortner_full.xy 1.xy
17 cp cmsky_cortner_full.rdls 1.rd
18 cp cmsky_cortner_full.wcs 1.wcs
19 cp cmsky_cortner_full.jpg 1.jpg
20 # or:
21 # wget "http://live.astrometry.net/status.php?job=alpha-201003-01883980&get=wcs.fits" -O 1.wcs
22 # wget "http://live.astrometry.net/status.php?job=alpha-201003-01883980&get=index.rd.fits" -O 1.rd
23 verify -w 1.wcs -x 1.xy -r 1.rd -v
24
25 **/
26
27 #include "starutil.h"
28 #include "mathutil.h"
29 #include "bl.h"
30 #include "matchobj.h"
31 #include "xylist.h"
32 #include "rdlist.h"
33 #include "ioutils.h"
34 #include "starkd.h"
35 #include "boilerplate.h"
36 #include "sip.h"
37 #include "sip_qfits.h"
38 #include "log.h"
39 #include "fitsioutils.h"
40 #include "fit-wcs.h"
41 #include "verify.h"
42 #include "histogram2d.h"
43
44 #include "plotstuff.h"
45 #include "plotimage.h"
46 #include "cairoutils.h"
47
48 static const char* OPTIONS = "hx:w:r:vj:p:";
49
print_help(char * progname)50 void print_help(char* progname) {
51 BOILERPLATE_HELP_HEADER(stdout);
52 printf("\nUsage: %s\n"
53 " -w <WCS input file>\n"
54 " -x <xyls input file>\n"
55 " -r <rdls input file>\n"
56 " [-p <plot output file>]\n"
57 " [-v]: verbose\n"
58 " [-j <pixel-jitter>]: set pixel jitter (default 1.0)\n"
59 "\n", progname);
60 }
61
62
main(int argc,char ** args)63 int main(int argc, char** args) {
64 int c;
65 char* xylsfn = NULL;
66 char* wcsfn = NULL;
67 char* rdlsfn = NULL;
68 char* plotfn = NULL;
69
70 xylist_t* xyls = NULL;
71 rdlist_t* rdls = NULL;
72 sip_t sip;
73 int i;
74 int W, H;
75 double pixeljitter = 1.0;
76 int loglvl = LOG_MSG;
77 double wcsscale;
78
79 fits_use_error_system();
80
81 while ((c = getopt(argc, args, OPTIONS)) != -1) {
82 switch (c) {
83 case 'p':
84 plotfn = optarg;
85 break;
86 case 'j':
87 pixeljitter = atof(optarg);
88 break;
89 case 'h':
90 print_help(args[0]);
91 exit(0);
92 case 'r':
93 rdlsfn = optarg;
94 break;
95 case 'x':
96 xylsfn = optarg;
97 break;
98 case 'w':
99 wcsfn = optarg;
100 break;
101 case 'v':
102 loglvl++;
103 break;
104 }
105 }
106 if (optind != argc) {
107 print_help(args[0]);
108 exit(-1);
109 }
110 if (!xylsfn || !wcsfn || !rdlsfn) {
111 print_help(args[0]);
112 exit(-1);
113 }
114 log_init(loglvl);
115
116 // read WCS.
117 logmsg("Trying to parse SIP header from %s...\n", wcsfn);
118 if (!sip_read_header_file(wcsfn, &sip)) {
119 logmsg("Failed to parse SIP header from %s.\n", wcsfn);
120 }
121 // image W, H
122 W = sip.wcstan.imagew;
123 H = sip.wcstan.imageh;
124 if ((W == 0.0) || (H == 0.0)) {
125 logmsg("WCS file %s didn't contain IMAGEW and IMAGEH headers.\n", wcsfn);
126 // FIXME - use bounds of xylist?
127 exit(-1);
128 }
129 wcsscale = sip_pixel_scale(&sip);
130 logmsg("WCS scale: %g arcsec/pixel\n", wcsscale);
131
132 // read XYLS.
133 xyls = xylist_open(xylsfn);
134 if (!xyls) {
135 logmsg("Failed to read an xylist from file %s.\n", xylsfn);
136 exit(-1);
137 }
138
139 // read RDLS.
140 rdls = rdlist_open(rdlsfn);
141 if (!rdls) {
142 logmsg("Failed to read an rdlist from file %s.\n", rdlsfn);
143 exit(-1);
144 }
145
146 {
147 // (x,y) positions of field stars.
148 double* fieldpix;
149 int Nfield;
150 double* indexpix;
151 starxy_t* xy;
152 rd_t* rd;
153 int Nindex;
154
155 xy = xylist_read_field(xyls, NULL);
156 if (!xy) {
157 logmsg("Failed to read xyls entries.\n");
158 exit(-1);
159 }
160 Nfield = starxy_n(xy);
161 fieldpix = starxy_to_xy_array(xy, NULL);
162 logmsg("Found %i field objects\n", Nfield);
163
164 // Project RDLS into pixel space.
165 rd = rdlist_read_field(rdls, NULL);
166 if (!rd) {
167 logmsg("Failed to read rdls entries.\n");
168 exit(-1);
169 }
170 Nindex = rd_n(rd);
171 logmsg("Found %i indx objects\n", Nindex);
172 indexpix = malloc(2 * Nindex * sizeof(double));
173 for (i=0; i<Nindex; i++) {
174 anbool ok;
175 double ra = rd_getra(rd, i);
176 double dec = rd_getdec(rd, i);
177 ok = sip_radec2pixelxy(&sip, ra, dec, indexpix + i*2, indexpix + i*2 + 1);
178 assert(ok);
179 }
180
181 logmsg("CRPIX is (%g,%g)\n", sip.wcstan.crpix[0], sip.wcstan.crpix[1]);
182
183 {
184 double* fieldsigma2s = malloc(Nfield * sizeof(double));
185 int besti;
186 int* theta;
187 double logodds;
188 double Q2, R2;
189 double qc[2];
190 double gamma;
191
192 // HACK -- quad radius-squared
193 Q2 = square(100.0);
194 qc[0] = sip.wcstan.crpix[0];
195 qc[1] = sip.wcstan.crpix[1];
196 // HACK -- variance growth rate wrt radius.
197 gamma = 1.0;
198
199 for (i=0; i<Nfield; i++) {
200 R2 = distsq(qc, fieldpix + 2*i, 2);
201 fieldsigma2s[i] = square(pixeljitter) * (1.0 + gamma * R2/Q2);
202 }
203
204 logodds = verify_star_lists(indexpix, Nindex,
205 fieldpix, fieldsigma2s, Nfield,
206 W*H,
207 0.25,
208 log(1e-100),
209 log(1e100),
210 &besti, NULL, &theta, NULL);
211
212 logmsg("Logodds: %g\n", logodds);
213
214 if (TRUE) {
215 for (i=0; i<Nfield; i++) {
216 if (theta[i] < 0)
217 continue;
218 printf("%g %g %g %g\n", fieldpix[2*i+0], fieldpix[2*i+1],
219 rd_getra(rd, theta[i]), rd_getdec(rd, theta[i]));
220 }
221 }
222
223 if (plotfn) {
224 plot_args_t pargs;
225 plotimage_t* img;
226 cairo_t* cairo;
227
228 plotstuff_init(&pargs);
229 pargs.outformat = PLOTSTUFF_FORMAT_PNG;
230 pargs.outfn = plotfn;
231 img = plotstuff_get_config(&pargs, "image");
232 img->format = PLOTSTUFF_FORMAT_JPG;
233 plot_image_set_filename(img, "1.jpg");
234 plot_image_setsize(&pargs, img);
235 plotstuff_run_command(&pargs, "image");
236 cairo = pargs.cairo;
237 // red circles around every field star.
238 cairo_set_color(cairo, "red");
239 for (i=0; i<Nfield; i++) {
240 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_CIRCLE,
241 fieldpix[2*i+0], fieldpix[2*i+1],
242 2.0 * sqrt(fieldsigma2s[i]));
243 cairo_stroke(cairo);
244 }
245 // green crosshairs at every index star.
246 cairo_set_color(cairo, "green");
247 for (i=0; i<Nindex; i++) {
248 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_XCROSSHAIR,
249 indexpix[2*i+0], indexpix[2*i+1],
250 3);
251 cairo_stroke(cairo);
252 }
253
254 // thick white circles for corresponding field stars.
255 cairo_set_line_width(cairo, 2);
256 for (i=0; i<Nfield; i++) {
257 if (theta[i] < 0)
258 continue;
259 cairo_set_color(cairo, "white");
260 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_CIRCLE,
261 fieldpix[2*i+0], fieldpix[2*i+1],
262 2.0 * sqrt(fieldsigma2s[i]));
263 cairo_stroke(cairo);
264 // thick cyan crosshairs for corresponding index stars.
265 cairo_set_color(cairo, "cyan");
266 cairoutils_draw_marker(cairo, CAIROUTIL_MARKER_XCROSSHAIR,
267 indexpix[2*theta[i]+0],
268 indexpix[2*theta[i]+1],
269 3);
270 cairo_stroke(cairo);
271
272 }
273
274 plotstuff_output(&pargs);
275 }
276
277 free(theta);
278 free(fieldsigma2s);
279 }
280
281 free(fieldpix);
282 free(indexpix);
283 }
284
285
286
287 if (xylist_close(xyls)) {
288 logmsg("Failed to close XYLS file.\n");
289 }
290 return 0;
291 }
292
