1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5
6 /**
7 * Accepts an image or xylist, plus command-line options, and produces
8 * an augmented xylist.
9 */
10
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <string.h>
14 #include <ctype.h>
15 #include <sys/types.h>
16 #include <sys/wait.h>
17 #include <sys/types.h>
18 #include <sys/stat.h>
19 #include <unistd.h>
20 #include <libgen.h>
21 #include <getopt.h>
22 #include <math.h>
23 #include <assert.h>
24
25 #include "os-features.h"
26 #include "ioutils.h"
27 #include "fileutils.h"
28 #include "bl.h"
29 #include "an-bool.h"
30 #include "solver.h"
31 #include "fitsioutils.h"
32 #include "solverutils.h"
33 #include "sip_qfits.h"
34 #include "tabsort.h"
35 #include "cut-table.h"
36 #include "errors.h"
37 #include "fits-guess-scale.h"
38 #include "image2xy-files.h"
39 #include "resort-xylist.h"
40 #include "an-opts.h"
41 #include "augment-xylist.h"
42 #include "log.h"
43 #include "anqfits.h"
44
45 static void delete_existing_an_headers(qfits_header* hdr);
46
augment_xylist_init(augment_xylist_t * axy)47 void augment_xylist_init(augment_xylist_t* axy) {
48 memset(axy, 0, sizeof(augment_xylist_t));
49 axy->tempdir = NULL;
50 axy->tweak = TRUE;
51 axy->tweakorder = 2;
52 axy->depths = il_new(4);
53 axy->fields = il_new(16);
54 axy->verifywcs = sl_new(4);
55 axy->verifywcs_ext = il_new(4);
56 axy->tagalong = sl_new(4);
57 axy->try_verify = TRUE;
58 axy->resort = TRUE;
59 axy->ra_center = HUGE_VAL;
60 axy->dec_center = HUGE_VAL;
61 axy->parity = PARITY_BOTH;
62 axy->uniformize = 10;
63 axy->verify_uniformize = TRUE;
64 }
65
augment_xylist_free_contents(augment_xylist_t * axy)66 void augment_xylist_free_contents(augment_xylist_t* axy) {
67 sl_free2(axy->verifywcs);
68 il_free(axy->verifywcs_ext);
69 sl_free2(axy->tagalong);
70 il_free(axy->depths);
71 il_free(axy->fields);
72 if (axy->predistort)
73 sip_free(axy->predistort);
74 }
75
augment_xylist_print_special_opts(an_option_t * opt,bl * opts,int index,FILE * fid,void * extra)76 void augment_xylist_print_special_opts(an_option_t* opt, bl* opts, int index,
77 FILE* fid, void* extra) {
78 if (!strcmp(opt->name, "image")) {
79 fprintf(fid, "%s",
80 " ( -i / --image <image-input-file>\n"
81 " OR -x / --xylist <xylist-input-file> ): input file\n");
82 } else if (!strcmp(opt->name, "scale-units")) {
83 fprintf(fid, "%s",
84 " -u / --scale-units <units>: in what units are the lower and upper bounds?\n"
85 " choices: \"degwidth\", \"degw\", \"dw\" : width of the image, in degrees (default)\n"
86 " \"arcminwidth\", \"amw\", \"aw\" : width of the image, in arcminutes\n"
87 " \"arcsecperpix\", \"app\": arcseconds per pixel\n"
88 " \"focalmm\": 35-mm (width-based) equivalent focal length\n"
89 );
90 } else if (!strcmp(opt->name, "depth")) {
91 fprintf(fid, "%s",
92 " -d / --depth <number or range>: number of field objects to look at, or range\n"
93 " of numbers; 1 is the brightest star, so \"-d 10\" or \"-d 1-10\" mean look\n"
94 " at the top ten brightest stars only.\n");
95 } else if (!strcmp(opt->name, "xylist-only")) {
96 fprintf(fid, "%s",
97 "The following options are valid for xylist inputs only:\n");
98 } else if (!strcmp(opt->name, "fields")) {
99 fprintf(fid, "%s",
100 " -F / --fields <number or range>: the FITS extension(s) to solve, inclusive\n");
101 } else if (!strcmp(opt->name, "options")) {
102 fprintf(fid, "Options include:\n");
103 }
104 }
105
106 static an_option_t options[] = {
107 {'i', "image", required_argument, NULL, NULL},
108 {'x', "xylist", required_argument, NULL, NULL},
109 {'o', "out", required_argument, "filename",
110 "output augmented xylist filename"},
111 {'\x01', "options", no_argument, NULL, NULL},
112 {'h', "help", no_argument, NULL,
113 "print this help message" },
114 {'v', "verbose", no_argument, NULL,
115 "be more chatty" },
116 {'7', "no-delete-temp", no_argument, NULL,
117 "don't delete temp files (for debugging)\n"},
118 {'L', "scale-low", required_argument, "scale",
119 "lower bound of image scale estimate"},
120 {'H', "scale-high", required_argument, "scale",
121 "upper bound of image scale estimate"},
122 {'u', "scale-units", required_argument, "units", NULL},
123 {'8', "parity", required_argument, "pos/neg",
124 "only check for matches with positive/negative parity (default: try both)"},
125 {'c', "code-tolerance", required_argument, "distance",
126 "matching distance for quads (default: 0.01)"},
127 {'E', "pixel-error", required_argument, "pixels",
128 "for verification, size of pixel positional error (default: 1)"},
129 {'q', "quad-size-min", required_argument, "fraction",
130 "minimum size of quads to try, as a fraction of the smaller image dimension, default: 0.1"},
131 {'Q', "quad-size-max", required_argument, "fraction",
132 "maximum size of quads to try, as a fraction of the image hypotenuse, default 1.0"},
133 {'\x82', "odds-to-tune-up", required_argument, "odds",
134 "odds ratio at which to try tuning up a match that isn't good enough to solve (default: 1e6)"},
135 {'[', "odds-to-solve", required_argument, "odds",
136 "odds ratio at which to consider a field solved (default: 1e9)"},
137 {'#', "odds-to-reject", required_argument, "odds",
138 "odds ratio at which to reject a hypothesis (default: 1e-100)"},
139 {'%', "odds-to-stop-looking", required_argument, "odds",
140 "odds ratio at which to stop adding stars when evaluating a hypothesis (default: HUGE_VAL)"},
141 {'^', "use-source-extractor", no_argument, NULL,
142 "use SourceExtractor rather than built-in image2xy to find sources"},
143 {'&', "source-extractor-config", required_argument, "filename",
144 "use the given SourceExtractor config file. "
145 "Note that CATALOG_NAME and CATALOG_TYPE values will be over-ridden by command-line values. "
146 "This option implies --use-source-extractor."},
147 {'*', "source-extractor-path", required_argument, "filename",
148 "use the given path to the SourceExtractor executable. Default: just 'source-extractor', assumed to be in your PATH."
149 " Note that you can give command-line args here too (but put them in quotes), eg: --source-extractor-path 'source-extractor -DETECT_TYPE CCD'. "
150 "This option implies --use-source-extractor."},
151 {'3', "ra", required_argument, "degrees or hh:mm:ss",
152 "only search in indexes within 'radius' of the field center given by 'ra' and 'dec'"},
153 {'4', "dec", required_argument, "degrees or [+-]dd:mm:ss",
154 "only search in indexes within 'radius' of the field center given by 'ra' and 'dec'"},
155 {'5', "radius", required_argument, "degrees",
156 "only search in indexes within 'radius' of the field center given by ('ra', 'dec')"},
157 {'d', "depth", required_argument, NULL, NULL},
158 {'|', "objs", required_argument, "int",
159 "cut the source list to have this many items (after sorting, if applicable)."},
160 {'l', "cpulimit", required_argument, "seconds",
161 "give up solving after the specified number of seconds of CPU time"},
162 {'r', "resort", no_argument, NULL,
163 "sort the star brightnesses by background-subtracted flux; the default is to sort using a"
164 "compromise between background-subtracted and non-background-subtracted flux"},
165 {'6', "extension", required_argument, "int",
166 "FITS extension to read image from."},
167 {';', "invert", no_argument, NULL,
168 "invert the image (for black-on-white images)"},
169 {'z', "downsample", required_argument, "int",
170 "downsample the image by factor <int> before running source extraction"},
171 {']', "no-background-subtraction", no_argument, NULL,
172 "don't try to estimate a smoothly-varying sky background during source extraction."},
173 {'{', "sigma", required_argument, "float",
174 "set the noise level in the image"},
175 {'\x93', "nsigma", required_argument, "float",
176 "number of sigma for a source detection; default 8"},
177 {'9', "no-remove-lines", no_argument, NULL,
178 "don't remove horizontal and vertical overdensities of sources."},
179 {':', "uniformize", required_argument, "int",
180 "select sources uniformly using roughly this many boxes (0=disable; default 10)"},
181 {'\x81', "no-verify-uniformize", no_argument, NULL,
182 "don't uniformize the field stars during verification"},
183 {'\x83', "no-verify-dedup", no_argument, NULL,
184 "don't deduplicate the field stars during verification"},
185 {'C', "cancel", required_argument, "filename",
186 "filename whose creation signals the process to stop"},
187 {'S', "solved", required_argument, "filename",
188 "output file to mark that the solver succeeded"},
189 {'I', "solved-in", required_argument, "filename",
190 "input filename for solved file"},
191 {'M', "match", required_argument, "filename",
192 "output filename for match file"},
193 {'R', "rdls", required_argument, "filename",
194 "output filename for RDLS file"},
195 {'\x80', "sort-rdls", required_argument, "column",
196 "sort the RDLS file by this column; default is ascending; use "
197 "\"-column\" to sort \"column\" in descending order instead."},
198 {'}',"tag", required_argument, "column",
199 "grab tag-along column from index into RDLS file"},
200 {'<', "tag-all", no_argument, NULL,
201 "grab all tag-along columns from index into RDLS file"},
202 {'j', "scamp-ref", required_argument, "filename",
203 "output filename for SCAMP reference catalog"},
204 {'B', "corr", required_argument, "filename",
205 "output filename for correspondences"},
206 {'W', "wcs", required_argument, "filename",
207 "output filename for WCS file"},
208 {'P', "pnm", required_argument, "filename",
209 "save the PNM file as <filename>"},
210 {'f', "force-ppm", no_argument, NULL,
211 "force the PNM file to be a PPM"},
212 {'k', "keep-xylist", required_argument, "filename",
213 "save the (unaugmented) xylist to <filename>"},
214 {'A', "dont-augment", no_argument, NULL,
215 "quit after writing the unaugmented xylist"},
216 {'V', "verify", required_argument, "filename",
217 "try to verify an existing WCS file"},
218 {'\x92', "verify-ext", required_argument, "extension",
219 "HDU from which to read WCS to verify; set this BEFORE --verify"},
220 {'y', "no-verify", no_argument, NULL,
221 "ignore existing WCS headers in FITS input images"},
222 {'g', "guess-scale", no_argument, NULL,
223 "try to guess the image scale from the FITS headers"},
224 {'>', "crpix-center", no_argument, NULL,
225 "set the WCS reference point to the image center"},
226 {'/', "crpix-x", required_argument, "pix",
227 "set the WCS reference point to the given position"},
228 {'\\', "crpix-y", required_argument, "pix",
229 "set the WCS reference point to the given position"},
230 {'T', "no-tweak", no_argument, NULL,
231 "don't fine-tune WCS by computing a SIP polynomial"},
232 {'t', "tweak-order", required_argument, "int",
233 "polynomial order of SIP WCS corrections"},
234 {'\x86', "predistort", required_argument, "filename",
235 "apply the inverse distortion in this SIP WCS header before solving"},
236 {'\x94', "xscale", required_argument, "factor",
237 "for rectangular pixels: factor to apply to measured X positions to make pixels square"},
238 {'m', "temp-dir", required_argument, "dir",
239 "where to put temp files, default /tmp"},
240 // placeholder for printing "The following are for xylist inputs only"
241 {'\0', "xylist-only", no_argument, NULL, NULL},
242 {'F', "fields", required_argument, NULL, NULL},
243 {'w', "width", required_argument, "pixels",
244 "specify the field width"},
245 {'e', "height", required_argument, "pixels",
246 "specify the field height"},
247 {'X', "x-column", required_argument, "column-name",
248 "the FITS column containing the X coordinate of the sources"},
249 {'Y', "y-column", required_argument, "column-name",
250 "the FITS column containing the Y coordinate of the sources"},
251 {'s', "sort-column", required_argument, "column-name",
252 "the FITS column that should be used to sort the sources"},
253 {'a', "sort-ascending", no_argument, NULL,
254 "sort in ascending order (smallest first); default is descending order"},
255 };
256
augment_xylist_print_help(FILE * fid)257 void augment_xylist_print_help(FILE* fid) {
258 bl* opts;
259 opts = opts_from_array(options, sizeof(options)/sizeof(an_option_t), NULL);
260 opts_print_help(opts, fid, augment_xylist_print_special_opts, NULL);
261 bl_free(opts);
262 }
263
augment_xylist_add_options(bl * opts)264 void augment_xylist_add_options(bl* opts) {
265 bl* myopts = opts_from_array(options, sizeof(options)/sizeof(an_option_t), NULL);
266 bl_append_list(opts, myopts);
267 bl_free(myopts);
268 }
269
270 static int parse_fields_string(il* fields, const char* str);
271
augment_xylist_parse_option(char argchar,char * optarg,augment_xylist_t * axy)272 int augment_xylist_parse_option(char argchar, char* optarg,
273 augment_xylist_t* axy) {
274 double d;
275 int verify_extension = -1;
276 //printf("parsing option %c (%i)\n", argchar, (int)argchar);
277 switch (argchar) {
278 case '\x80':
279 axy->sort_rdls = optarg;
280 break;
281 case '\x81':
282 axy->verify_uniformize = FALSE;
283 break;
284 case '\x82':
285 axy->odds_to_tune_up = atof(optarg);
286 break;
287 case '\x83':
288 axy->verify_dedup = FALSE;
289 break;
290 case '\x86':
291 axy->predistort = sip_read_header_file(optarg, NULL);
292 if (!axy->predistort) {
293 ERROR("Failed to read SIP header file \"%s\" for pre-distortion values", optarg);
294 return -1;
295 }
296 break;
297 case '\x94':
298 axy->pixel_xscale = atof(optarg);
299 break;
300 case ';':
301 axy->invert_image = TRUE;
302 break;
303 case '>':
304 axy->set_crpix_center = TRUE;
305 axy->set_crpix = TRUE;
306 break;
307 case '/':
308 axy->crpix[0] = atof(optarg);
309 axy->set_crpix = TRUE;
310 break;
311 case '\\':
312 axy->crpix[1] = atof(optarg);
313 axy->set_crpix = TRUE;
314 break;
315 case '}':
316 sl_append(axy->tagalong, optarg);
317 break;
318 case '<':
319 axy->tagalong_all = TRUE;
320 break;
321 case '{':
322 axy->image_sigma = atof(optarg);
323 break;
324 case '\x93':
325 axy->image_nsigma = atof(optarg);
326 break;
327 case '^':
328 axy->use_source_extractor = TRUE;
329 break;
330 case '&':
331 axy->source_extractor_config = optarg;
332 axy->use_source_extractor = TRUE;
333 break;
334 case '*':
335 axy->source_extractor_path = optarg;
336 axy->use_source_extractor = TRUE;
337 break;
338 case '9':
339 axy->no_removelines = TRUE;
340 break;
341 case ':':
342 axy->uniformize = atoi(optarg);
343 break;
344 case '3':
345 axy->ra_center = atora(optarg);
346 if (axy->ra_center == HUGE_VAL) {
347 ERROR("Couldn't understand your RA center argument \"%s\"", optarg);
348 return -1;
349 }
350 break;
351 case '4':
352 axy->dec_center = atodec(optarg);
353 if (axy->dec_center == HUGE_VAL) {
354 ERROR("Couldn't understand your Dec center argument \"%s\"", optarg);
355 return -1;
356 }
357 break;
358 case '5':
359 axy->search_radius = atof(optarg);
360 break;
361 case '6':
362 axy->extension = atoi(optarg);
363 break;
364 case '[':
365 axy->odds_to_solve = atof(optarg);
366 break;
367 case '#':
368 axy->odds_to_bail = atof(optarg);
369 break;
370 case '%':
371 axy->odds_to_stoplooking = atof(optarg);
372 break;
373 case ']':
374 axy->no_bg_subtraction = TRUE;
375 break;
376 case '8':
377 if (streq(optarg, "pos")) {
378 axy->parity = PARITY_NORMAL;
379 } else if (streq(optarg, "neg")) {
380 axy->parity = PARITY_FLIP;
381 } else {
382 ERROR("Couldn't understand your Parity argument \"%s\": must be \"pos\" or \"neg\"", optarg);
383 return -1;
384 }
385 break;
386 case 'B':
387 axy->corrfn = optarg;
388 break;
389 case 'y':
390 axy->try_verify = FALSE;
391 break;
392 case 'q':
393 d = atof(optarg);
394 if (d < 0.0 || d > 1.0) {
395 ERROR("quad size fraction (-q / --quad-size-min) must be between 0.0 and 1.0");
396 return -1;
397 }
398 axy->quadsize_min = d;
399 break;
400 case 'Q':
401 d = atof(optarg);
402 if (d < 0.0 || d > 1.0) {
403 ERROR("quad size fraction (-Q / --quad-size-max) must be between 0.0 and 1.0");
404 return -1;
405 }
406 axy->quadsize_max = d;
407 break;
408 case 'l':
409 axy->cpulimit = atof(optarg);
410 break;
411 case 'A':
412 axy->dont_augment = TRUE;
413 break;
414 case 'z':
415 axy->downsample = atoi(optarg);
416 break;
417 case 'r':
418 axy->resort = FALSE;
419 break;
420 case 'E':
421 axy->pixelerr = atof(optarg);
422 break;
423 case 'c':
424 axy->codetol = atof(optarg);
425 break;
426 case 'v':
427 axy->verbosity++;
428 break;
429 case '7':
430 axy->no_delete_temp = TRUE;
431 break;
432 case '\x92':
433 verify_extension = atoi(optarg);
434 break;
435 case 'V':
436 sl_append(axy->verifywcs, optarg);
437 il_append(axy->verifywcs_ext,
438 (verify_extension >= 0 ? verify_extension : axy->extension));
439 break;
440 case 'I':
441 axy->solvedinfn = optarg;
442 break;
443 case 'k':
444 axy->keepxylsfn = optarg;
445 break;
446 case 's':
447 axy->sortcol = optarg;
448 axy->resort = FALSE;
449 break;
450 case 'a':
451 axy->sort_ascending = TRUE;
452 break;
453 case 'X':
454 axy->xcol = optarg;
455 break;
456 case 'Y':
457 axy->ycol = optarg;
458 break;
459 case 'm':
460 axy->tempdir = optarg;
461 break;
462 case 'F':
463 if (parse_fields_string(axy->fields, optarg)) {
464 ERROR("Failed to parse fields specification \"%s\"", optarg);
465 return -1;
466 }
467 break;
468 case 'd':
469 if (parse_depth_string(axy->depths, optarg)) {
470 ERROR("Failed to parse depth specification: \"%s\"", optarg);
471 return -1;
472 }
473 break;
474 case '|':
475 axy->cutobjs = atoi(optarg);
476 break;
477 case 'o':
478 axy->axyfn = optarg;
479 break;
480 case 'i':
481 axy->imagefn = optarg;
482 break;
483 case 'x':
484 axy->xylsfn = optarg;
485 break;
486 case 'L':
487 axy->scalelo = atof(optarg);
488 break;
489 case 'H':
490 axy->scalehi = atof(optarg);
491 break;
492 case 'u':
493 axy->scaleunit = parse_scale_units(optarg);
494 if (axy->scaleunit == -1) {
495 ERROR("Unknown image scale units \"%s\". See \"solve-field -h\" for allowed values.\n", optarg);
496 return -1;
497 }
498 break;
499 case 'w':
500 axy->W = atoi(optarg);
501 break;
502 case 'e':
503 axy->H = atoi(optarg);
504 break;
505 case 'T':
506 axy->tweak = FALSE;
507 break;
508 case 't':
509 axy->tweakorder = atoi(optarg);
510 break;
511 case 'P':
512 axy->pnmfn = optarg;
513 break;
514 case 'f':
515 axy->force_ppm = TRUE;
516 break;
517 case 'g':
518 axy->guess_scale = TRUE;
519 break;
520 case 'S':
521 axy->solvedfn = optarg;
522 break;
523 case 'C':
524 axy->cancelfn = optarg;
525 break;
526 case 'M':
527 axy->matchfn = optarg;
528 break;
529 case 'R':
530 axy->rdlsfn = optarg;
531 break;
532 case 'j':
533 axy->scampfn = optarg;
534 break;
535 case 'W':
536 axy->wcsfn = optarg;
537 break;
538 default:
539 return 1;
540 }
541
542 return 0;
543 }
544
parse_scale_units(const char * units)545 int parse_scale_units(const char* units) {
546 if (strcaseeq(units, "degwidth") ||
547 strcaseeq(units, "degw") ||
548 strcaseeq(units, "dw")) {
549 return SCALE_UNITS_DEG_WIDTH;
550 } else if (strcaseeq(units, "arcminwidth") ||
551 strcaseeq(units, "amw") ||
552 strcaseeq(units, "aw")) {
553 return SCALE_UNITS_ARCMIN_WIDTH;
554 } else if (strcaseeq(units, "arcsecperpix") ||
555 strcaseeq(units, "app")) {
556 return SCALE_UNITS_ARCSEC_PER_PIX;
557 } else if (strcaseeq(units, "focalmm")) {
558 return SCALE_UNITS_FOCAL_MM;
559 }
560 return -1;
561 }
562
563 // run(): ppmtopgm, pnmtofits, source_extractor
564 // backtick(): pnmfile, image2pnm
565
append_escape(sl * list,const char * fn)566 static void append_escape(sl* list, const char* fn) {
567 sl_append_nocopy(list, shell_escape(fn));
568 }
append_executable(sl * list,const char * fn,const char * me)569 static void append_executable(sl* list, const char* fn, const char* me) {
570 char* exec = find_executable(fn, me);
571 if (!exec) {
572 char* binfn = NULL;
573 asprintf_safe(&binfn, "../bin/%s", fn);
574 exec = find_executable(binfn, me);
575 free(binfn);
576 }
577 if (!exec) {
578 ERROR("Couldn't find executable \"%s\"", fn);
579 exit(-1);
580 }
581 sl_append_nocopy(list, shell_escape(exec));
582 free(exec);
583 }
584
backtick(sl * cmd,anbool verbose)585 static sl* backtick(sl* cmd, anbool verbose) {
586 char* cmdstr = sl_implode(cmd, " ");
587 sl* lines;
588 logverb("Running: %s\n", cmdstr);
589 if (run_command_get_outputs(cmdstr, &lines, NULL)) {
590 ERROR("Failed to run command: %s", cmdstr);
591 free(cmdstr);
592 exit(-1);
593 }
594 free(cmdstr);
595 sl_remove_all(cmd);
596 return lines;
597 }
598
run(sl * cmd,anbool verbose)599 static void run(sl* cmd, anbool verbose) {
600 if (verbose) {
601 char* cmdstr = sl_implode(cmd, " ");
602 logverb("Running: %s\n", cmdstr);
603 if (run_command_get_outputs(cmdstr, NULL, NULL)) {
604 ERROR("Failed to run command: %s", cmdstr);
605 free(cmdstr);
606 exit(-1);
607 }
608 free(cmdstr);
609 sl_remove_all(cmd);
610 } else {
611 sl* lines = backtick(cmd, verbose);
612 if (verbose) {
613 int i;
614 for (i=0; i<sl_size(lines); i++)
615 logverb(" %s\n", sl_get(lines, i));
616 }
617 sl_free2(lines);
618 }
619 }
620
add_sip_coeffs(qfits_header * hdr,const char * prefix,const sip_t * sip)621 static void add_sip_coeffs(qfits_header* hdr, const char* prefix, const sip_t* sip) {
622 char key[64];
623 int m, n, order;
624
625 if (sip->a_order) {
626 sprintf(key, "%sSAO", prefix);
627 order = sip->a_order;
628 fits_header_add_int(hdr, key, order, "SIP forward polynomial order");
629 for (m=0; m<=order; m++) {
630 for (n=0; (m+n)<=order; n++) {
631 if (m+n < 1)
632 continue;
633 sprintf(key, "%sA%i%i", prefix, m, n);
634 fits_header_add_double(hdr, key, sip->a[m][n], "");
635 sprintf(key, "%sB%i%i", prefix, m, n);
636 fits_header_add_double(hdr, key, sip->b[m][n], "");
637 }
638 }
639 }
640 if (sip->ap_order) {
641 order = sip->ap_order;
642 sprintf(key, "%sSAPO", prefix);
643 fits_header_add_int(hdr, key, order, "SIP reverse polynomial order");
644 for (m=0; m<=order; m++) {
645 for (n=0; (m+n)<=order; n++) {
646 if (m+n < 1)
647 continue;
648 sprintf(key, "%sAP%i%i", prefix, m, n);
649 fits_header_add_double(hdr, key, sip->ap[m][n], "");
650 sprintf(key, "%sBP%i%i", prefix, m, n);
651 fits_header_add_double(hdr, key, sip->bp[m][n], "");
652 }
653 }
654 }
655 }
656
augment_xylist(augment_xylist_t * axy,const char * me)657 int augment_xylist(augment_xylist_t* axy,
658 const char* me) {
659 // tempfiles to delete when we finish
660 sl* tempfiles;
661 sl* cmd;
662 anbool verbose = axy->verbosity > 0;
663 int i, I;
664 //anbool guessed_scale = FALSE;
665 anbool dosort = FALSE;
666 char* xylsfn;
667 qfits_header* hdr = NULL;
668 anbool addwh = TRUE;
669 FILE* fout = NULL;
670 char *fitsimgfn = NULL;
671 dl* scales;
672 char* nolinesfn = NULL;
673 char* sortedxylsfn = NULL;
674 char* unixylsfn = NULL;
675 char* cutxylsfn = NULL;
676 // as we process the image (uncompress it, eg), keep track of extension
677 // (along with axy->fitsimgfn if keep_fitsimg is set)
678 axy->fitsimgext = axy->extension;
679
680 cmd = sl_new(16);
681 tempfiles = sl_new(4);
682 scales = dl_new(4);
683
684 if (axy->imagefn) {
685 // if --image is given:
686 // -run image2pnm
687 // -if it's a FITS image, keep the original
688 // -otherwise, run ppmtopgm (if necessary) and pnmtofits.
689 // -run image2xy to generate xylist
690 char *uncompressedfn;
691 char *pnmfn = NULL;
692 sl* lines;
693 anbool iscompressed = FALSE;
694 char* line;
695 char pnmtype;
696 char typestr[256];
697 anbool want_pnm = TRUE;
698
699 uncompressedfn = create_temp_file("uncompressed", axy->tempdir);
700 sl_append_nocopy(tempfiles, uncompressedfn);
701
702 if (axy->assume_fits_image) {
703 axy->isfits = TRUE;
704 if (!axy->pnmfn) {
705 qfits_header* hdr;
706 want_pnm = FALSE;
707 // We need to get image W,H from the FITS header.
708 logverb("Reading FITS image \"%s\" to find image size\n", axy->imagefn);
709 hdr = anqfits_get_header2(axy->imagefn, axy->extension);
710 axy->W = qfits_header_getint(hdr, "NAXIS1", -1);
711 axy->H = qfits_header_getint(hdr, "NAXIS2", -1);
712 qfits_header_destroy(hdr);
713 if (axy->W == -1 || axy->H == -1) {
714 ERROR("Failed to find size of FITS image \"%s\": got NAXIS1 = %i, NAXIS2 = %i\n",
715 axy->imagefn, axy->W, axy->H);
716 return -1;
717 }
718 logverb(" got FITS image size %i x %i\n", axy->W, axy->H);
719 }
720 }
721
722 if (want_pnm) {
723 if (axy->pnmfn)
724 pnmfn = axy->pnmfn;
725 else {
726 pnmfn = create_temp_file("pnm", axy->tempdir);
727 sl_append_nocopy(tempfiles, pnmfn);
728 }
729
730 append_executable(cmd, "image2pnm", me);
731 if (axy->extension) {
732 sl_append(cmd, "--extension");
733 sl_appendf(cmd, "%i", axy->extension);
734 }
735 sl_append(cmd, "--infile");
736 append_escape(cmd, axy->imagefn);
737 sl_append(cmd, "--uncompressed-outfile");
738 append_escape(cmd, uncompressedfn);
739 sl_append(cmd, "--outfile");
740 append_escape(cmd, pnmfn);
741 if (axy->force_ppm)
742 sl_append(cmd, "--ppm");
743 sl_append(cmd, "--mydir");
744 append_escape(cmd, me);
745 lines = backtick(cmd, verbose);
746 axy->isfits = FALSE;
747 for (i=0; i<sl_size(lines); i++) {
748 logverb(" %s\n", sl_get(lines, i));
749 if (streq("fits", sl_get(lines, i)))
750 axy->isfits = TRUE;
751 if (streq("compressed", sl_get(lines, i)))
752 iscompressed = TRUE;
753 if (streq("fz", sl_get(lines, i)))
754 axy->fitsimgext = 0;
755 }
756 sl_free2(lines);
757
758 // Get image W, H, depth.
759 sl_append(cmd, "pnmfile");
760 append_escape(cmd, pnmfn);
761
762 lines = backtick(cmd, verbose);
763
764 if (sl_size(lines) == 0) {
765 ERROR("Got no output from the \"pnmfile\" program.");
766 exit(-1);
767 }
768 line = sl_get(lines, 0);
769 // eg "/tmp/pnm: PPM raw, 800 by 510 maxval 255"
770 if (strlen(pnmfn) + 1 >= strlen(line)) {
771 ERROR("Failed to parse output from pnmfile: \"%s\"", line);
772 exit(-1);
773 }
774 line += strlen(pnmfn) + 1;
775 // "PBM raw, 750 by 864"
776 if (sscanf(line, " P%cM %255s %d by %d",
777 &pnmtype, typestr, &axy->W, &axy->H) != 4) {
778 ERROR("Failed to parse output from pnmfile: \"%s\"\n", line);
779 exit(-1);
780 }
781 sl_free2(lines);
782 }
783
784 if (axy->isfits) {
785
786 if (iscompressed)
787 fitsimgfn = uncompressedfn;
788 else
789 fitsimgfn = axy->imagefn;
790
791 if (axy->try_verify) {
792 char* errstr;
793 sip_t sip;
794 anbool ok;
795 // Try to read WCS header from FITS image; if successful,
796 // add it to the list of WCS headers to verify.
797 logverb("Looking for a WCS header in FITS input image %s ext %i\n", fitsimgfn, axy->fitsimgext);
798
799 // FIXME - Right now we just try to read SIP/TAN -
800 // obviously this should be more flexible and robust.
801 errors_start_logging_to_string();
802 memset(&sip, 0, sizeof(sip_t));
803 ok = (sip_read_header_file_ext(fitsimgfn, axy->fitsimgext, &sip) != NULL);
804 errstr = errors_stop_logging_to_string(": ");
805 if (ok) {
806 logmsg("Found an existing WCS header, will try to verify it.\n");
807 sl_append(axy->verifywcs, fitsimgfn);
808 il_append(axy->verifywcs_ext, axy->fitsimgext);
809 } else {
810 logverb("Failed to read a SIP or TAN header from FITS image.\n");
811 logverb(" (reason: %s)\n", errstr);
812 }
813 free(errstr);
814 }
815
816 } else {
817 fitsimgfn = create_temp_file("fits", axy->tempdir);
818 sl_append_nocopy(tempfiles, fitsimgfn);
819
820 if (pnmtype == 'P') {
821 logverb("Converting PPM image to FITS...\n");
822
823 sl_append(cmd, "ppmtopgm");
824 append_escape(cmd, pnmfn);
825 sl_append(cmd, "|");
826 append_executable(cmd, "an-pnmtofits", me);
827 sl_append(cmd, ">");
828 append_escape(cmd, fitsimgfn);
829
830 run(cmd, verbose);
831
832 } else if (pnmtype == 'G') {
833 logverb("Converting PGM image to FITS...\n");
834
835 append_executable(cmd, "an-pnmtofits", me);
836 append_escape(cmd, pnmfn);
837 sl_append(cmd, ">");
838 append_escape(cmd, fitsimgfn);
839
840 run(cmd, verbose);
841
842 } else if (pnmtype == 'B') {
843 logverb("Converting PBM image to FITS...\n");
844
845 sl_append(cmd, "pbmtopgm 3 3");
846 append_escape(cmd, pnmfn);
847 sl_append(cmd, "|");
848 append_executable(cmd, "an-pnmtofits", me);
849 sl_append(cmd, ">");
850 append_escape(cmd, fitsimgfn);
851 run(cmd, verbose);
852 } else {
853 assert(0);
854 }
855
856 }
857
858 if (axy->keep_fitsimg) {
859 axy->fitsimgfn = strdup(fitsimgfn);
860 sl_remove_string(tempfiles, fitsimgfn);
861 }
862
863 logmsg("Extracting sources...\n");
864 xylsfn = create_temp_file("xyls", axy->tempdir);
865 sl_append_nocopy(tempfiles, xylsfn);
866
867 if (axy->use_source_extractor) {
868 if (axy->source_extractor_path)
869 sl_append(cmd, axy->source_extractor_path);
870 else
871 sl_append(cmd, "source-extractor");
872
873 if (axy->source_extractor_config)
874 sl_appendf(cmd, "-c %s", axy->source_extractor_config);
875 else {
876 char* paramfn;
877 char* paramstr;
878 char* filterfn;
879 char* filterstr;
880
881 paramfn = create_temp_file("param", axy->tempdir);
882 sl_append_nocopy(tempfiles, paramfn);
883 paramstr = "X_IMAGE\nY_IMAGE\nMAG_AUTO\nFLUX_AUTO";
884 if (write_file(paramfn, paramstr, strlen(paramstr))) {
885 ERROR("Failed to write Source Extractor parameters to temp file \"%s\"", paramfn);
886 exit(-1);
887 }
888 sl_appendf(cmd, "-PARAMETERS_NAME %s", paramfn);
889 if (verbose)
890 sl_append(cmd, "-VERBOSE_TYPE FULL");
891
892 axy->xcol = "X_IMAGE";
893 axy->ycol = "Y_IMAGE";
894 axy->sortcol = "MAG_AUTO";
895 axy->sort_ascending = TRUE;
896
897 filterfn = create_temp_file("filter", axy->tempdir);
898 sl_append_nocopy(tempfiles, filterfn);
899 filterstr = "CONV NORM\n"
900 "# 5x5 convolution mask of a gaussian PSF with FWHM = 2.0 pixels.\n"
901 "0.006319 0.040599 0.075183 0.040599 0.006319\n"
902 "0.040599 0.260856 0.483068 0.260856 0.040599\n"
903 "0.075183 0.483068 0.894573 0.483068 0.075183\n"
904 "0.040599 0.260856 0.483068 0.260856 0.040599\n"
905 "0.006319 0.040599 0.075183 0.040599 0.006319\n";
906 if (write_file(filterfn, filterstr, strlen(filterstr))) {
907 ERROR("Failed to write Source Extractor convolution filter to temp file \"%s\"", filterfn);
908 exit(-1);
909 }
910 sl_appendf(cmd, "-FILTER_NAME %s", filterfn);
911 }
912
913 sl_append(cmd, "-CATALOG_TYPE FITS_1.0");
914 sl_appendf(cmd, "-CATALOG_NAME %s", xylsfn);
915 append_escape(cmd, fitsimgfn);
916
917 logverb("Running Source Extractor: output file is %s\n", xylsfn);
918 run(cmd, verbose);
919
920 } else {
921 simplexy_t sxyparams;
922 logverb("Running image2xy: input=%s, output=%s, ext=%i\n", fitsimgfn, xylsfn, axy->fitsimgext);
923
924 // we have to delete the temp file because otherwise image2xy is too timid to overwrite it.
925 if (unlink(xylsfn)) {
926 SYSERROR("Failed to delete temp file %s", xylsfn);
927 exit(-1);
928 }
929
930 memset(&sxyparams, 0, sizeof(simplexy_t));
931 // The other params get set to defaults for float or u8 images.
932 sxyparams.nobgsub = axy->no_bg_subtraction;
933 sxyparams.sigma = axy->image_sigma;
934 sxyparams.invert = axy->invert_image;
935 sxyparams.plim = axy->image_nsigma;
936
937 // MAGIC 3: downsample by a factor of 2, up to 3 times.
938 if (image2xy_files(fitsimgfn, xylsfn, TRUE, axy->downsample, 3,
939 axy->fitsimgext, 0, &sxyparams)) {
940 ERROR("Source extraction failed");
941 exit(-1);
942 }
943 }
944 dosort = TRUE;
945
946 } else {
947 // xylist.
948 xylsfn = axy->xylsfn;
949 if (axy->sortcol)
950 dosort = TRUE;
951
952 if (axy->extension && (axy->extension != 1)) {
953 // Copy just this extension to a temp file.
954 FILE* fout;
955 FILE* fin;
956 off_t offset;
957 off_t nbytes;
958 off_t nprimary;
959 anqfits_t* anq;
960 char* extfn;
961
962 anq = anqfits_open_hdu(xylsfn, axy->extension);
963 if (!anq) {
964 ERROR("Failed to open xyls file %s up to extension %i",
965 xylsfn, axy->extension);
966 exit(-1);
967 }
968
969 fin = fopen(xylsfn, "rb");
970 if (!fin) {
971 SYSERROR("Failed to open xyls file \"%s\"", xylsfn);
972 exit(-1);
973 }
974
975 nprimary = anqfits_header_size(anq, 0) + anqfits_data_size(anq, 0);
976 offset = anqfits_header_start(anq, axy->extension);
977 nbytes = anqfits_header_size(anq, axy->extension) +
978 anqfits_data_size(anq, axy->extension);
979 anqfits_close(anq);
980 extfn = create_temp_file("ext", axy->tempdir);
981 sl_append_nocopy(tempfiles, extfn);
982 fout = fopen(extfn, "wb");
983 if (!fout) {
984 SYSERROR("Failed to open temp file \"%s\" to write extension",
985 extfn);
986 exit(-1);
987 }
988 logverb("Copying ext %i of %s to temp %s\n", axy->extension,
989 xylsfn, extfn);
990 if (pipe_file_offset(fin, 0, nprimary, fout)) {
991 ERROR("Failed to copy the primary HDU of xylist file %s to %s",
992 xylsfn, extfn);
993 exit(-1);
994 }
995 if (pipe_file_offset(fin, offset, nbytes, fout)) {
996 ERROR("Failed to copy HDU %i of xylist file %s to %s",
997 axy->extension, xylsfn, extfn);
998 exit(-1);
999 }
1000 fclose(fin);
1001 if (fclose(fout)) {
1002 SYSERROR("Failed to close %s", extfn);
1003 exit(-1);
1004 }
1005 xylsfn = extfn;
1006 }
1007 }
1008
1009 if (axy->guess_scale && (fitsimgfn || !axy->imagefn)) {
1010 dl* estscales = NULL;
1011 char* infn = (fitsimgfn ? fitsimgfn : xylsfn);
1012 fits_guess_scale(infn, NULL, &estscales);
1013 for (i=0; i<dl_size(estscales); i++) {
1014 double scale = dl_get(estscales, i);
1015 logverb("Scale estimate: %g\n", scale);
1016 dl_append(scales, scale * 0.99);
1017 dl_append(scales, scale * 1.01);
1018 //guessed_scale = TRUE;
1019 }
1020 dl_free(estscales);
1021 }
1022
1023 // remove lines
1024 // sort
1025 // uniformize
1026 // cut
1027 if (axy->keepxylsfn) {
1028 // Figure out which is the last stage to run, and set its output
1029 // file to "keepxylsfn".
1030 if (axy->cutobjs) {
1031 cutxylsfn = axy->keepxylsfn;
1032 } else if (axy->uniformize) {
1033 unixylsfn = axy->keepxylsfn;
1034 } else if (dosort) {
1035 sortedxylsfn = axy->keepxylsfn;
1036 } else if (!axy->no_removelines) {
1037 nolinesfn = axy->keepxylsfn;
1038 } else {
1039 // copy xylsfn to axy->keepxylsfn.
1040 if (copy_file(xylsfn, axy->keepxylsfn)) {
1041 ERROR("Failed to copy xyls file \"%s\" to \"%s\"",
1042 xylsfn, axy->keepxylsfn);
1043 return -1;
1044 }
1045 }
1046 }
1047
1048 if (!axy->no_removelines) {
1049 if (!nolinesfn) {
1050 nolinesfn = create_temp_file("removelines", axy->tempdir);
1051 sl_append_nocopy(tempfiles, nolinesfn);
1052 }
1053 logverb("Removing lines of (spurious) sources from xylist \"%s\", writing to \"%s\"\n",
1054 xylsfn, nolinesfn);
1055 append_executable(cmd, "removelines", me);
1056 if (axy->xcol)
1057 sl_appendf(cmd, "-X %s", axy->xcol);
1058 if (axy->ycol)
1059 sl_appendf(cmd, "-Y %s", axy->ycol);
1060 //if (axy->extension)
1061 // sl_appendf(cmd, "-e %i", axy->extension);
1062 append_escape(cmd, xylsfn);
1063 append_escape(cmd, nolinesfn);
1064 run(cmd, verbose);
1065 xylsfn = nolinesfn;
1066 }
1067
1068 if (dosort) {
1069 anbool do_tabsort = FALSE;
1070
1071 if (!axy->sortcol)
1072 axy->sortcol = "FLUX";
1073 if (!axy->bgcol)
1074 axy->bgcol = "BACKGROUND";
1075
1076 if (!sortedxylsfn) {
1077 sortedxylsfn = create_temp_file("sorted", axy->tempdir);
1078 sl_append_nocopy(tempfiles, sortedxylsfn);
1079 }
1080
1081 if (axy->resort) {
1082 char* err;
1083 int rtn;
1084 logverb("Sorting file \"%s\" to \"%s\" using columns flux (%s) and background (%s), %sscending\n",
1085 xylsfn, sortedxylsfn, axy->sortcol, axy->bgcol, axy->sort_ascending?"a":"de");
1086 errors_start_logging_to_string();
1087 rtn = resort_xylist(xylsfn, sortedxylsfn, axy->sortcol, axy->bgcol, axy->sort_ascending);
1088 err = errors_stop_logging_to_string(": ");
1089 if (rtn) {
1090 logmsg("Sorting brightness using %s and BACKGROUND columns failed; falling back to %s.\n",
1091 axy->sortcol, axy->sortcol);
1092 logverb("Reason: %s\n", err);
1093 do_tabsort = TRUE;
1094 }
1095 free(err);
1096
1097 } else
1098 do_tabsort = TRUE;
1099
1100 if (do_tabsort) {
1101 logverb("Sorting by brightness: input=%s, output=%s, column=%s.\n",
1102 xylsfn, sortedxylsfn, axy->sortcol);
1103 tabsort(xylsfn, sortedxylsfn,
1104 axy->sortcol, !axy->sort_ascending);
1105 }
1106 xylsfn = sortedxylsfn;
1107 }
1108
1109 if (axy->uniformize) {
1110 if (!unixylsfn) {
1111 unixylsfn = create_temp_file("uniform", axy->tempdir);
1112 sl_append_nocopy(tempfiles, unixylsfn);
1113 }
1114 append_executable(cmd, "uniformize", me);
1115 sl_appendf(cmd, "-n %i", axy->uniformize);
1116 if (axy->xcol)
1117 sl_appendf(cmd, "-X %s", axy->xcol);
1118 if (axy->ycol)
1119 sl_appendf(cmd, "-Y %s", axy->ycol);
1120 //if (axy->extension)
1121 // sl_appendf(cmd, "-e %i", axy->extension);
1122 append_escape(cmd, xylsfn);
1123 append_escape(cmd, unixylsfn);
1124 run(cmd, verbose);
1125 xylsfn = unixylsfn;
1126 }
1127
1128 if (axy->cutobjs) {
1129 // cut the source lists to at most "cutobjs" objects.
1130 if (!cutxylsfn) {
1131 cutxylsfn = create_temp_file("cut", axy->tempdir);
1132 sl_append_nocopy(tempfiles, cutxylsfn);
1133 }
1134
1135 if (cut_table(xylsfn, cutxylsfn, axy->cutobjs)) {
1136 ERROR("Failed to cut table %s to %i entries; output file %s", xylsfn, axy->cutobjs, cutxylsfn);
1137 return -1;
1138 }
1139 xylsfn = cutxylsfn;
1140 }
1141
1142 if (axy->dont_augment)
1143 // done!
1144 goto cleanup;
1145
1146 // start piling FITS headers in there.
1147 hdr = anqfits_get_header2(xylsfn, 0);
1148 if (!hdr) {
1149 ERROR("Failed to read FITS header from file %s", xylsfn);
1150 exit(-1);
1151 }
1152
1153 // delete any existing processing directives
1154 delete_existing_an_headers(hdr);
1155
1156 if (!(axy->W && axy->H)) {
1157 // Look for existing IMAGEW and IMAGEH in primary header.
1158 axy->W = qfits_header_getint(hdr, "IMAGEW", 0);
1159 axy->H = qfits_header_getint(hdr, "IMAGEH", 0);
1160 if (axy->W && axy->H) {
1161 addwh = FALSE;
1162 } else {
1163 // Look for IMAGEW and IMAGEH headers in first extension, else bail.
1164 qfits_header* hdr2 = anqfits_get_header2(xylsfn, 1);
1165 axy->W = qfits_header_getint(hdr2, "IMAGEW", 0);
1166 axy->H = qfits_header_getint(hdr2, "IMAGEH", 0);
1167 qfits_header_destroy(hdr2);
1168 }
1169 if (!(axy->W && axy->H)) {
1170 ERROR("Error: image width and height must be specified for XYLS inputs");
1171 exit(-1);
1172 }
1173 }
1174
1175 // we may write long filenames.
1176 fits_header_add_longstring_boilerplate(hdr);
1177
1178 if (addwh) {
1179 fits_header_add_int(hdr, "IMAGEW", axy->W, "image width");
1180 fits_header_add_int(hdr, "IMAGEH", axy->H, "image height");
1181 }
1182 qfits_header_add(hdr, "ANRUN", "T", "Solve this field!", NULL);
1183
1184 if (axy->cpulimit > 0)
1185 fits_header_add_double(hdr, "ANCLIM", axy->cpulimit, "CPU time limit (seconds)");
1186
1187 if (axy->xcol)
1188 qfits_header_add(hdr, "ANXCOL", axy->xcol, "Name of column containing X coords", NULL);
1189 if (axy->ycol)
1190 qfits_header_add(hdr, "ANYCOL", axy->ycol, "Name of column containing Y coords", NULL);
1191
1192 if (axy->tagalong_all)
1193 qfits_header_add(hdr, "ANTAGALL", "T", "Tag-along all columns from index to RDLS", NULL);
1194 else
1195 for (i=0; i<sl_size(axy->tagalong); i++) {
1196 char key[64];
1197 sprintf(key, "ANTAG%i", i+1);
1198 qfits_header_add(hdr, key, sl_get(axy->tagalong, i), "Tag-along column from index to RDLS", NULL);
1199 }
1200
1201 if (axy->sort_rdls)
1202 qfits_header_add(hdr, "ANRDSORT", axy->sort_rdls, "Sort RDLS file by this column", NULL);
1203
1204 qfits_header_add(hdr, "ANVERUNI", axy->verify_uniformize ? "T":"F", "Uniformize field during verification", NULL);
1205 qfits_header_add(hdr, "ANVERDUP", axy->verify_dedup ? "T":"F", "Deduplicate field during verification", NULL);
1206
1207 if (axy->odds_to_tune_up)
1208 fits_header_add_double(hdr, "ANODDSTU", axy->odds_to_tune_up, "Odds ratio to tune up a match");
1209 if (axy->odds_to_solve)
1210 fits_header_add_double(hdr, "ANODDSSL", axy->odds_to_solve, "Odds ratio to consider a field solved");
1211 if (axy->odds_to_bail)
1212 fits_header_add_double(hdr, "ANODDSBL", axy->odds_to_bail, "Odds ratio to consider a hypothesis rejected");
1213 if (axy->odds_to_stoplooking)
1214 fits_header_add_double(hdr, "ANODDSST", axy->odds_to_stoplooking, "Odds ratio to stop trying to improve the odds ratio");
1215
1216 if ((axy->scalelo > 0.0) || (axy->scalehi > 0.0)) {
1217 double appu, appl;
1218 switch (axy->scaleunit) {
1219 case SCALE_UNITS_DEG_WIDTH:
1220 logverb("Scale range: %g to %g degrees wide\n", axy->scalelo, axy->scalehi);
1221 appl = deg2arcsec(axy->scalelo) / (double)axy->W;
1222 appu = deg2arcsec(axy->scalehi) / (double)axy->W;
1223 logverb("Image width %i pixels; arcsec per pixel range %g %g\n", axy->W, appl, appu);
1224 break;
1225 case SCALE_UNITS_ARCMIN_WIDTH:
1226 logverb("Scale range: %g to %g arcmin wide\n", axy->scalelo, axy->scalehi);
1227 appl = arcmin2arcsec(axy->scalelo) / (double)axy->W;
1228 appu = arcmin2arcsec(axy->scalehi) / (double)axy->W;
1229 logverb("Image width %i pixels; arcsec per pixel range %g %g\n", axy->W, appl, appu);
1230 break;
1231 case SCALE_UNITS_ARCSEC_PER_PIX:
1232 logverb("Scale range: %g to %g arcsec/pixel\n", axy->scalelo, axy->scalehi);
1233 appl = axy->scalelo;
1234 appu = axy->scalehi;
1235 break;
1236 case SCALE_UNITS_FOCAL_MM:
1237 logverb("Scale range: %g to %g mm focal length\n", axy->scalelo, axy->scalehi);
1238 // "35 mm" film is 36 mm wide.
1239 appu = rad2arcsec(atan(36. / (2. * axy->scalelo))) / (double)axy->W;
1240 appl = rad2arcsec(atan(36. / (2. * axy->scalehi))) / (double)axy->W;
1241 logverb("Image width %i pixels; arcsec per pixel range %g %g\n", axy->W, appl, appu);
1242 break;
1243 default:
1244 ERROR("Unknown scale unit code %i\n", axy->scaleunit);
1245 return -1;
1246 }
1247 dl_append(scales, appl);
1248 dl_append(scales, appu);
1249 }
1250
1251 /* Hmm, do we want this??
1252 if ((dl_size(axy->scales) > 0) && guessed_scale)
1253 qfits_header_add(hdr, "ANAPPDEF", "T", "try the default scale range too.", NULL);
1254 */
1255
1256 for (i=0; i<dl_size(scales)/2; i++) {
1257 char key[64];
1258 double lo = dl_get(scales, 2*i);
1259 double hi = dl_get(scales, 2*i + 1);
1260 if (lo > 0.0) {
1261 sprintf(key, "ANAPPL%i", i+1);
1262 fits_header_add_double(hdr, key, lo, "scale: arcsec/pixel min");
1263 }
1264 if (hi > 0.0) {
1265 sprintf(key, "ANAPPU%i", i+1);
1266 fits_header_add_double(hdr, key, hi, "scale: arcsec/pixel max");
1267 }
1268 }
1269
1270 if (axy->quadsize_min > 0.0)
1271 fits_header_add_double(hdr, "ANQSFMIN", axy->quadsize_min, "minimum quad size: fraction");
1272 if (axy->quadsize_max > 0.0)
1273 fits_header_add_double(hdr, "ANQSFMAX", axy->quadsize_max, "maximum quad size: fraction");
1274
1275 if (axy->set_crpix) {
1276 if (axy->set_crpix_center) {
1277 qfits_header_add(hdr, "ANCRPIXC", "T", "Set CRPIX to the image center.", NULL);
1278 } else {
1279 fits_header_add_double(hdr, "ANCRPIX1", axy->crpix[0], "Set CRPIX1 to this val.");
1280 fits_header_add_double(hdr, "ANCRPIX2", axy->crpix[1], "Set CRPIX2 to this val.");
1281 }
1282 }
1283
1284 qfits_header_add(hdr, "ANTWEAK", (axy->tweak ? "T" : "F"), (axy->tweak ? "Tweak: yes please!" : "Tweak: no, thanks."), NULL);
1285 if (axy->tweak && axy->tweakorder)
1286 fits_header_add_int(hdr, "ANTWEAKO", axy->tweakorder, "Tweak order");
1287
1288 if (axy->solvedfn)
1289 fits_header_addf_longstring(hdr, "ANSOLVED", "solved output file", "%s", axy->solvedfn);
1290 if (axy->solvedinfn)
1291 fits_header_addf_longstring(hdr, "ANSOLVIN", "solved input file", "%s", axy->solvedinfn);
1292 if (axy->cancelfn)
1293 fits_header_addf_longstring(hdr, "ANCANCEL", "cancel output file", "%s", axy->cancelfn);
1294 if (axy->matchfn)
1295 fits_header_addf_longstring(hdr, "ANMATCH", "match output file", "%s", axy->matchfn);
1296 if (axy->rdlsfn)
1297 fits_header_addf_longstring(hdr, "ANRDLS", "ra-dec output file", "%s", axy->rdlsfn);
1298 if (axy->scampfn)
1299 fits_header_addf_longstring(hdr, "ANSCAMP", "SCAMP reference catalog output file", "%s", axy->scampfn);
1300 if (axy->wcsfn)
1301 fits_header_addf_longstring(hdr, "ANWCS", "WCS header output filename", "%s", axy->wcsfn);
1302 if (axy->corrfn)
1303 fits_header_addf_longstring(hdr, "ANCORR", "Correspondences output filename", "%s", axy->corrfn);
1304 if (axy->codetol > 0.0)
1305 fits_header_add_double(hdr, "ANCTOL", axy->codetol, "code tolerance");
1306 if (axy->pixelerr > 0.0)
1307 fits_header_add_double(hdr, "ANPOSERR", axy->pixelerr, "star pos'n error (pixels)");
1308
1309 if (axy->parity != PARITY_BOTH) {
1310 if (axy->parity == PARITY_NORMAL)
1311 qfits_header_add(hdr, "ANPARITY", "POS", "det(CD) > 0", NULL);
1312 else if (axy->parity == PARITY_FLIP)
1313 qfits_header_add(hdr, "ANPARITY", "NEG", "det(CD) < 0", NULL);
1314 }
1315
1316 if ((axy->ra_center != HUGE_VAL) &&
1317 (axy->dec_center != HUGE_VAL) &&
1318 (axy->search_radius >= 0.0)) {
1319 fits_header_add_double(hdr, "ANERA", axy->ra_center, "RA center estimate (deg)");
1320 fits_header_add_double(hdr, "ANEDEC", axy->dec_center, "Dec center estimate (deg)");
1321 fits_header_add_double(hdr, "ANERAD", axy->search_radius, "Search radius from estimated posn (deg)");
1322 }
1323
1324 for (i=0; i<il_size(axy->depths)/2; i++) {
1325 int depthlo, depthhi;
1326 char key[64];
1327 depthlo = il_get(axy->depths, 2*i);
1328 depthhi = il_get(axy->depths, 2*i + 1);
1329 sprintf(key, "ANDPL%i", (i+1));
1330 fits_header_addf(hdr, key, "", "%i", depthlo);
1331 sprintf(key, "ANDPU%i", (i+1));
1332 fits_header_addf(hdr, key, "", "%i", depthhi);
1333 }
1334
1335 for (i=0; i<il_size(axy->fields)/2; i++) {
1336 int lo = il_get(axy->fields, 2*i);
1337 int hi = il_get(axy->fields, 2*i + 1);
1338 char key[64];
1339 if (lo == hi) {
1340 sprintf(key, "ANFD%i", (i+1));
1341 fits_header_add_int(hdr, key, lo, "field to solve");
1342 } else {
1343 sprintf(key, "ANFDL%i", (i+1));
1344 fits_header_add_int(hdr, key, lo, "field range: low");
1345 sprintf(key, "ANFDU%i", (i+1));
1346 fits_header_add_int(hdr, key, hi, "field range: high");
1347 }
1348 }
1349
1350 I = 0;
1351 for (i=0; i<sl_size(axy->verifywcs); i++) {
1352 sip_t sip;
1353 const char* fn;
1354 int ext;
1355
1356 fn = sl_get(axy->verifywcs, i);
1357 ext = il_get(axy->verifywcs_ext, i);
1358 if (!sip_read_header_file_ext(fn, ext, &sip)) {
1359 ERROR("Failed to parse WCS header from file \"%s\" ext %i", fn, ext);
1360 continue;
1361 }
1362
1363 I++;
1364 {
1365 tan_t* wcs = &(sip.wcstan);
1366 // note, this initialization has to happen *after* you read the WCS header :)
1367 double vals[] = { wcs->crval[0], wcs->crval[1],
1368 wcs->crpix[0], wcs->crpix[1],
1369 wcs->cd[0][0], wcs->cd[0][1],
1370 wcs->cd[1][0], wcs->cd[1][1] };
1371 char key[64];
1372 char* keys[] = { "ANW%iPIX1", "ANW%iPIX2", "ANW%iVAL1", "ANW%iVAL2",
1373 "ANW%iCD11", "ANW%iCD12", "ANW%iCD21", "ANW%iCD22" };
1374 int j;
1375 for (j = 0; j < 8; j++) {
1376 sprintf(key, keys[j], I);
1377 fits_header_add_double(hdr, key, vals[j], "");
1378 }
1379
1380 sprintf(key, "ANW%i", I);
1381 add_sip_coeffs(hdr, key, &sip);
1382 }
1383 }
1384
1385 if (axy->predistort) {
1386 fits_header_add_double(hdr, "ANDPIX0", axy->predistort->wcstan.crpix[0], "Pre-distortion ref pix x");
1387 fits_header_add_double(hdr, "ANDPIX1", axy->predistort->wcstan.crpix[1], "Pre-distortion ref pix y");
1388 add_sip_coeffs(hdr, "AND", axy->predistort);
1389 }
1390
1391 if (axy->pixel_xscale) {
1392 fits_header_add_double(hdr, "ANPXSCAL", axy->pixel_xscale, "x scaling to make square pixels");
1393 }
1394
1395 fout = fopen(axy->axyfn, "wb");
1396 if (!fout) {
1397 SYSERROR("Failed to open output file %s", axy->axyfn);
1398 exit(-1);
1399 }
1400
1401 logverb("Writing headers to file %s\n", axy->axyfn);
1402
1403 if (qfits_header_dump(hdr, fout)) {
1404 ERROR("Failed to write FITS header");
1405 exit(-1);
1406 }
1407 qfits_header_destroy(hdr);
1408
1409 // copy blocks from xyls to output.
1410 {
1411 FILE* fin;
1412 off_t offset;
1413 off_t nbytes;
1414 anqfits_t* anq;
1415 int ext = 1;
1416
1417 anq = anqfits_open_hdu(xylsfn, ext);
1418 if (!anq) {
1419 ERROR("Failed to open xyls file %s up to extension %i", xylsfn, ext);
1420 exit(-1);
1421 }
1422 offset = anqfits_header_start(anq, ext);
1423 nbytes = anqfits_header_size(anq, ext) + anqfits_data_size(anq, ext);
1424
1425 logverb("Copying data block of file %s to output %s.\n",
1426 xylsfn, axy->axyfn);
1427 anqfits_close(anq);
1428
1429 fin = fopen(xylsfn, "rb");
1430 if (!fin) {
1431 SYSERROR("Failed to open xyls file \"%s\"", xylsfn);
1432 exit(-1);
1433 }
1434
1435 if (pipe_file_offset(fin, offset, nbytes, fout)) {
1436 ERROR("Failed to copy the data segment of xylist file %s to %s",
1437 xylsfn, axy->axyfn);
1438 exit(-1);
1439 }
1440 fclose(fin);
1441 }
1442 fclose(fout);
1443
1444 cleanup:
1445 if (!axy->no_delete_temp) {
1446 for (i=0; i<sl_size(tempfiles); i++) {
1447 char* fn = sl_get(tempfiles, i);
1448 logverb("Deleting temp file %s\n", fn);
1449 if (unlink(fn)) {
1450 SYSERROR("Failed to delete temp file \"%s\"", fn);
1451 }
1452 }
1453 }
1454
1455 dl_free(scales);
1456 sl_free2(cmd);
1457 sl_free2(tempfiles);
1458 return 0;
1459 }
1460
delete_existing_an_headers(qfits_header * hdr)1461 static void delete_existing_an_headers(qfits_header* hdr) {
1462 int i,j,k;
1463 char key[64];
1464 qfits_header_del(hdr, "ANRUN");
1465 qfits_header_del(hdr, "ANCLIM");
1466 qfits_header_del(hdr, "ANXCOL");
1467 qfits_header_del(hdr, "ANYCOL");
1468 qfits_header_del(hdr, "ANTAGALL");
1469 for (i=0; i<100; i++) {
1470 sprintf(key, "ANTAG%i", i+1);
1471 qfits_header_del(hdr, key);
1472 sprintf(key, "ANAPPL%i", i+1);
1473 qfits_header_del(hdr, key);
1474 sprintf(key, "ANAPPU%i", i+1);
1475 qfits_header_del(hdr, key);
1476 sprintf(key, "ANDPL%i", i+1);
1477 qfits_header_del(hdr, key);
1478 sprintf(key, "ANDPU%i", i+1);
1479 qfits_header_del(hdr, key);
1480 sprintf(key, "ANFD%i", i+1);
1481 qfits_header_del(hdr, key);
1482 sprintf(key, "ANFDL%i", i+1);
1483 qfits_header_del(hdr, key);
1484 sprintf(key, "ANFDU%i", i+1);
1485 qfits_header_del(hdr, key);
1486 }
1487 for (i=0; i<10; i++) {
1488 sprintf(key, "ANW%iPIX1", i+1);
1489 qfits_header_del(hdr, key);
1490 sprintf(key, "ANW%iPIX2", i+1);
1491 qfits_header_del(hdr, key);
1492 sprintf(key, "ANW%iVAL1", i+1);
1493 qfits_header_del(hdr, key);
1494 sprintf(key, "ANW%iVAL2", i+1);
1495 qfits_header_del(hdr, key);
1496 sprintf(key, "ANW%iCD11", i+1);
1497 qfits_header_del(hdr, key);
1498 sprintf(key, "ANW%iCD12", i+1);
1499 qfits_header_del(hdr, key);
1500 sprintf(key, "ANW%iCD21", i+1);
1501 qfits_header_del(hdr, key);
1502 sprintf(key, "ANW%iCD22", i+1);
1503 qfits_header_del(hdr, key);
1504 sprintf(key, "ANW%iSAO", i+1);
1505 qfits_header_del(hdr, key);
1506 sprintf(key, "ANW%iSAPO", i+1);
1507 qfits_header_del(hdr, key);
1508 sprintf(key, "ANW%iSBO", i+1);
1509 qfits_header_del(hdr, key);
1510 sprintf(key, "ANW%iSBPO", i+1);
1511 qfits_header_del(hdr, key);
1512 for (j=0; j<10; j++) {
1513 for (k=0; k<10; k++) {
1514 sprintf(key, "ANW%iA%i%i", i+1, j, k);
1515 qfits_header_del(hdr, key);
1516 sprintf(key, "ANW%iB%i%i", i+1, j, k);
1517 qfits_header_del(hdr, key);
1518 sprintf(key, "ANW%iAP%i%i", i+1, j, k);
1519 qfits_header_del(hdr, key);
1520 sprintf(key, "ANW%iBP%i%i", i+1, j, k);
1521 qfits_header_del(hdr, key);
1522 }
1523 }
1524 }
1525 qfits_header_del(hdr, "ANRDSORT");
1526 qfits_header_del(hdr, "ANVERUNI");
1527 qfits_header_del(hdr, "ANVERDUP");
1528 qfits_header_del(hdr, "ANODDSTU");
1529 qfits_header_del(hdr, "ANODDSSL");
1530 qfits_header_del(hdr, "ANODDSBL");
1531 qfits_header_del(hdr, "ANODDSST");
1532 qfits_header_del(hdr, "ANQSFMIN");
1533 qfits_header_del(hdr, "ANQSFMAX");
1534 qfits_header_del(hdr, "ANCRPIXC");
1535 qfits_header_del(hdr, "ANCRPIX1");
1536 qfits_header_del(hdr, "ANCRPIX2");
1537 qfits_header_del(hdr, "ANTWEAK");
1538 qfits_header_del(hdr, "ANTWEAKO");
1539 qfits_header_del(hdr, "ANSOLVED");
1540 qfits_header_del(hdr, "ANSOLVIN");
1541 qfits_header_del(hdr, "ANCANCEL");
1542 qfits_header_del(hdr, "ANMATCH");
1543 qfits_header_del(hdr, "ANRDLS");
1544 qfits_header_del(hdr, "ANSCAMP");
1545 qfits_header_del(hdr, "ANWCS");
1546 qfits_header_del(hdr, "ANCORR");
1547 qfits_header_del(hdr, "ANCTOL");
1548 qfits_header_del(hdr, "ANPOSERR");
1549 qfits_header_del(hdr, "ANPARITY");
1550 qfits_header_del(hdr, "ANERA");
1551 qfits_header_del(hdr, "ANEDEC");
1552 qfits_header_del(hdr, "ANERAD");
1553 qfits_header_del(hdr, "ANDPIX0");
1554 qfits_header_del(hdr, "ANDPIX1");
1555 qfits_header_del(hdr, "ANDSAO");
1556 qfits_header_del(hdr, "ANDSBO");
1557 qfits_header_del(hdr, "ANDSAPO");
1558 qfits_header_del(hdr, "ANDSBPO");
1559 for (j=0; j<10; j++) {
1560 for (k=0; k<10; k++) {
1561 sprintf(key, "ANDA%i%i", j, k);
1562 qfits_header_del(hdr, key);
1563 sprintf(key, "ANDB%i%i", j, k);
1564 qfits_header_del(hdr, key);
1565 sprintf(key, "ANDAP%i%i", j, k);
1566 qfits_header_del(hdr, key);
1567 sprintf(key, "ANDBP%i%i", j, k);
1568 qfits_header_del(hdr, key);
1569 }
1570 }
1571 }
1572
parse_fields_string(il * fields,const char * str)1573 static int parse_fields_string(il* fields, const char* str) {
1574 // 10,11,20-25,30,40-50
1575 while (str && *str) {
1576 unsigned int lo, hi;
1577 int nread;
1578 if (sscanf(str, "%u-%u", &lo, &hi) == 2) {
1579 sscanf(str, "%*u-%*u%n", &nread);
1580 } else if (sscanf(str, "%u", &lo) == 1) {
1581 sscanf(str, "%*u%n", &nread);
1582 hi = lo;
1583 } else {
1584 ERROR("Failed to parse fragment: \"%s\"", str);
1585 return -1;
1586 }
1587 if (lo <= 0) {
1588 ERROR("Field number %i is invalid: must be >= 1.", lo);
1589 return -1;
1590 }
1591 if (lo > hi) {
1592 ERROR("Field range %i to %i is invalid: max must be >= min!", lo, hi);
1593 return -1;
1594 }
1595 il_append(fields, lo);
1596 il_append(fields, hi);
1597 str += nread;
1598 while ((*str == ',') || isspace((unsigned)(*str)))
1599 str++;
1600 }
1601 return 0;
1602 }
1603
1604