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 "os-features.h"
6 #include "starutil.h"
7 #include "mathutil.h"
8 #include "fit-wcs.h"
9 #include "xylist.h"
10 #include "rdlist.h"
11 #include "boilerplate.h"
12 #include "sip.h"
13 #include "sip_qfits.h"
14 #include "fitsioutils.h"
15 #include "anwcs.h"
16 #include "log.h"
17 #include "errors.h"
18
19 static const char* OPTIONS = "hw:e:tLx:y:W:H:N:o:v";
20
print_help(char * progname)21 void print_help(char* progname) {
22 BOILERPLATE_HELP_HEADER(stdout);
23 printf("\nUsage: %s\n"
24 " -w <WCS input file>\n"
25 " [-e <extension>] FITS HDU number to read WCS from (default 0 = primary)\n"
26 " [-t]: just use TAN projection, even if SIP extension exists.\n"
27 " [-L]: force WCSlib\n"
28 " [-x x-lo]\n"
29 " [-y y-lo]\n"
30 " [-W x-hi]\n"
31 " [-H y-hi]\n"
32 " [-N grid-n]\n"
33 " -o <WCS output file>\n"
34 " [-v]: verbose\n"
35 "\n", progname);
36 }
37
38
main(int argc,char ** args)39 int main(int argc, char** args) {
40 int c;
41 char* wcsfn = NULL;
42 char* outfn = NULL;
43 int ext = 0;
44 anbool forcetan = FALSE;
45 anbool forcewcslib = FALSE;
46 double xlo = 1;
47 double ylo = 1;
48 double xhi = HUGE_VAL;
49 double yhi = HUGE_VAL;
50 int N = 20;
51
52 double* xyz = NULL;
53 double* xy = NULL;
54
55 tan_t outwcs;
56 anwcs_t* inwcs = NULL;
57 int i, j;
58 double xstep, ystep;
59 int loglvl = LOG_MSG;
60
61 fits_use_error_system();
62
63 while ((c = getopt(argc, args, OPTIONS)) != -1) {
64 switch (c) {
65 case 'h':
66 print_help(args[0]);
67 exit(0);
68 case 'w':
69 wcsfn = optarg;
70 break;
71 case 'L':
72 forcewcslib = TRUE;
73 break;
74 case 't':
75 forcetan = TRUE;
76 break;
77 case 'o':
78 outfn = optarg;
79 break;
80 case 'e':
81 ext = atoi(optarg);
82 break;
83 case 'N':
84 N = atoi(optarg);
85 break;
86 case 'x':
87 xlo = atof(optarg);
88 break;
89 case 'y':
90 ylo = atof(optarg);
91 break;
92 case 'W':
93 xhi = atof(optarg);
94 break;
95 case 'H':
96 yhi = atof(optarg);
97 break;
98 case 'v':
99 loglvl++;
100 break;
101 }
102 }
103 if (optind != argc) {
104 print_help(args[0]);
105 exit(-1);
106 }
107 if (!wcsfn || !outfn) {
108 print_help(args[0]);
109 exit(-1);
110 }
111 log_init(loglvl);
112
113 // read WCS.
114 logmsg("Trying to read WCS header from \"%s\" ext %i...\n", wcsfn, ext);
115 if (forcewcslib) {
116 inwcs = anwcs_open_wcslib(wcsfn, ext);
117 } else if (forcetan) {
118 inwcs = anwcs_open_tan(wcsfn, ext);
119 } else {
120 inwcs = anwcs_open(wcsfn, ext);
121 }
122 if (!inwcs) {
123 ERROR("Failed to read WCS file \"%s\", extension %i", wcsfn, ext);
124 exit(-1);
125 }
126 logverb("Read WCS:\n");
127 if (log_get_level() >= LOG_VERB) {
128 anwcs_print(inwcs, log_get_fid());
129 }
130
131 if (xhi == HUGE_VAL) {
132 xhi = anwcs_imagew(inwcs);
133 logverb("Setting image width to %g\n", xhi);
134 }
135 if (yhi == HUGE_VAL) {
136 yhi = anwcs_imageh(inwcs);
137 logverb("Setting image height to %g\n", yhi);
138 }
139 // FIXME -- what if the user wants xhi or yhi == 0?
140 if (xhi == HUGE_VAL || xhi == 0) {
141 ERROR("Couldn't find the image size; please supply -W\n");
142 exit(-1);
143 }
144 if (yhi == HUGE_VAL || yhi == 0) {
145 ERROR("Couldn't find the image size; please supply -H\n");
146 exit(-1);
147 }
148
149 xstep = (xhi - xlo) / (N-1);
150 ystep = (yhi - ylo) / (N-1);
151
152 logverb("Evaluating WCS on a grid of %i x %i in X [%g,%g], Y [%g,%g]\n",
153 N, N, xlo, xhi, ylo, yhi);
154
155 xyz = (double*)malloc(sizeof(double) * 3 * N*N);
156 xy = (double*)malloc(sizeof(double) * 2 * N*N);
157 if (!xyz || !xy) {
158 ERROR("Failed to allocate %i xyz, xy coords", N*N);
159 exit(-1);
160 }
161
162 for (j=0; j<N; j++) {
163 for (i=0; i<N; i++) {
164 double x, y;
165 x = xlo + xstep * i;
166 y = ylo + ystep * j;
167 if (anwcs_pixelxy2xyz(inwcs, x, y, xyz + (j*N + i)*3)) {
168 ERROR("Failed to apply WCS to pixel coord (%g,%g)", x, y);
169 exit(-1);
170 }
171 xy[(j*N+i)*2 + 0] = x;
172 xy[(j*N+i)*2 + 1] = y;
173 }
174 }
175
176 logverb("Fitting TAN WCS\n");
177 if (fit_tan_wcs(xyz, xy, N*N, &outwcs, NULL)) {
178 ERROR("Failed to fit TAN WCS");
179 exit(-1);
180 }
181
182 if (tan_write_to_file(&outwcs, outfn)) {
183 ERROR("Failed to write TAN WCS header to file \"%s\"", outfn);
184 exit(-1);
185 }
186
187 free(xy);
188 free(xyz);
189 anwcs_free(inwcs);
190
191 return 0;
192 }
193