1 /*============================================================================
2 WCSLIB 7.7 - an implementation of the FITS WCS standard.
3 Copyright (C) 1995-2021, Mark Calabretta
4
5 This file is part of WCSLIB.
6
7 WCSLIB is free software: you can redistribute it and/or modify it under the
8 terms of the GNU Lesser General Public License as published by the Free
9 Software Foundation, either version 3 of the License, or (at your option)
10 any later version.
11
12 WCSLIB is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 more details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with WCSLIB. If not, see http://www.gnu.org/licenses.
19
20 Author: Mark Calabretta, Australia Telescope National Facility, CSIRO.
21 http://www.atnf.csiro.au/people/Mark.Calabretta
22 $Id: twcstab.c,v 7.7 2021/07/12 06:36:49 mcalabre Exp $
23 *=============================================================================
24 *
25 * twcstab tests wcstab() and also provides sample code for using it in
26 * conjunction with wcspih() and fits_read_wcstab(). Although this example
27 * and fits_read_wcstab() are based on the CFITSIO library, wcstab() itself
28 * is completely independent of it.
29 *
30 * The input FITS file is constructed by create_input() from a list of header
31 * keyrecords in wcstab.keyrec together with some hard-coded parameters. The
32 * output fits file, wcstab.fits, is left for inspection.
33 *
34 *===========================================================================*/
35
36 #include <math.h>
37 #include <stdio.h>
38 #include <string.h>
39
40 #include <fitsio.h>
41
42 #include <wcslib.h>
43 #include <getwcstab.h>
44
45 int create_input();
46 int do_wcs_stuff(fitsfile *fptr, struct wcsprm *wcs);
47 double lcprng();
48
main()49 int main()
50
51 {
52 char *header;
53 int i, nkeyrec, nreject, nwcs, stat[NWCSFIX], status = 0;
54 fitsfile *fptr;
55 struct wcsprm *wcs;
56
57
58 // Set line buffering in case stdout is redirected to a file, otherwise
59 // stdout and stderr messages will be jumbled (stderr is unbuffered).
60 setvbuf(stdout, NULL, _IOLBF, 0);
61
62 printf("Testing -TAB interpreter (twcstab.c)\n"
63 "------------------------------------\n\n");
64
65 // Create the input FITS test file.
66 if (create_input()) {
67 fprintf(stderr, "Failed to create FITS test file.");
68 return 1;
69 }
70
71 // Open the FITS test file and read the primary header.
72 fits_open_file(&fptr, "wcstab.fits", READONLY, &status);
73 if ((status = fits_hdr2str(fptr, 1, NULL, 0, &header, &nkeyrec,
74 &status))) {
75 fits_report_error(stderr, status);
76 return 1;
77 }
78
79
80 //-------------------------------------------------------------------------
81 // Basic steps required to interpret a FITS WCS header, including -TAB.
82 //-------------------------------------------------------------------------
83
84 // Parse the primary header of the FITS file.
85 if ((status = wcspih(header, nkeyrec, WCSHDR_all, 2, &nreject, &nwcs,
86 &wcs))) {
87 fprintf(stderr, "wcspih ERROR %d: %s.\n", status,wcshdr_errmsg[status]);
88 }
89
90 // Read coordinate arrays from the binary table extension.
91 if ((status = fits_read_wcstab(fptr, wcs->nwtb, (wtbarr *)wcs->wtb,
92 &status))) {
93 fits_report_error(stderr, status);
94 return 1;
95 }
96
97 // Translate non-standard WCS keyvalues.
98 if ((status = wcsfix(7, 0, wcs, stat))) {
99 for (i = 0; i < NWCSFIX; i++) {
100 if (stat[i] > 0) {
101 fprintf(stderr, "wcsfix ERROR %d: %s.\n", status,
102 wcsfix_errmsg[stat[i]]);
103 }
104 }
105
106 return 1;
107 }
108
109 //-------------------------------------------------------------------------
110 // The wcsprm struct is now ready for use.
111 //-------------------------------------------------------------------------
112
113 // Do something with it.
114 do_wcs_stuff(fptr, wcs);
115
116 // Finished with the FITS file.
117 fits_close_file(fptr, &status);
118 fits_free_memory(header, &status);
119
120 // Clean up.
121 status = wcsvfree(&nwcs, &wcs);
122
123 return 0;
124 }
125
126 //----------------------------------------------------------------------------
127
128 // The celestial plane is 256 x 128 with the table indexed at every fourth
129 // pixel on each axis, but the image is rotated by 5 deg so the table needs
130 // to be a bit wider than 65 x 33.
131
132 #define K1 70L
133 #define K2 40L
134
create_input()135 int create_input()
136
137 {
138 const double TWOPI = 2.0 * 3.14159265358979323846;
139
140 // These must match wcstab.keyrec.
141 const char infile[] = "test/wcstab.keyrec";
142 const long NAXIS1 = 256;
143 const long NAXIS2 = 128;
144 const long NAXIS3 = 4;
145 const char *ttype1[3] = {"CelCoords", "RAIndex", "DecIndex"};
146 const char *tform1[3] = {"5600E", "70E", "40E"};
147 const char *tunit1[3] = {"deg", "", ""};
148 const char *ttype2[4] = {"WaveIndex", "WaveCoord",
149 "TimeIndex", "TimeCoord"};
150 const char *tform2[4] = {"8E", "8D", "8E", "8D"};
151 const char *tunit2[4] = {"", "m", "", "a"};
152
153 // The remaining parameters may be chosen at will.
154 float refval[] = {150.0f, -30.0f};
155 float span[] = {5.0f, (5.0f*K2)/K1};
156 float sigma[] = {0.10f, 0.05f};
157 double wcoord[] = {0.21106114, 0.21076437, 2.0e-6, 2.2e-6,
158 500.0e-9, 650.0e-9, 1.24e-9, 2.48e-9};
159 double windex[] = {0.5, 1.5, 1.5, 2.5, 2.5, 3.5, 3.5, 4.5};
160 double tcoord[] = {1997.84512, 1997.84631, 1993.28451, 1993.28456,
161 2001.59234, 2001.59239, 2002.18265, 2002.18301};
162 double tindex[] = {0.0, 1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0};
163
164 char keyrec[84];
165 int status;
166 size_t iz;
167 long dummy, firstelem, k1, k2, p1, p2, p3;
168 float array[2*K1*K2], *fp, image[256];
169 double s, t, x1, x2, z, z1, z2;
170 FILE *stream;
171 fitsfile *fptr;
172
173 // Look for the input header keyrecords.
174 if ((stream = fopen(infile+5, "r")) == 0x0) {
175 if ((stream = fopen(infile, "r")) == 0x0) {
176 printf("ERROR opening %s\n", infile);
177 return 1;
178 }
179 }
180
181 // Create the FITS output file, deleting any pre-existing file.
182 status = 0;
183 fits_create_file(&fptr, "!wcstab.fits", &status);
184
185 // Convert header keyrecords to FITS.
186 while (fgets(keyrec, 82, stream) != NULL) {
187 // Ignore meta-comments (copyright information, etc.).
188 if (keyrec[0] == '#') continue;
189
190 // Strip off the newline.
191 iz = strlen(keyrec) - 1;
192 if (keyrec[iz] == '\n') keyrec[iz] = '\0';
193
194 fits_write_record(fptr, keyrec, &status);
195 }
196 fclose(stream);
197
198 // Create and write some phoney image data.
199 firstelem = 1;
200 for (p3 = 0; p3 < NAXIS3; p3++) {
201 for (p2 = 0; p2 < NAXIS2; p2++) {
202 fp = image;
203 s = (p3 + 1) * cos(0.2 * p2);
204 t = cos(0.8*p2);
205 for (p1 = 0; p1 < NAXIS1; p1++) {
206 // Do not adjust your set!
207 *(fp++) = sin(0.1*(p1+p2) + s) * cos(0.4*p1) * t;
208 }
209
210 fits_write_img_flt(fptr, 0L, firstelem, NAXIS1, image, &status);
211 firstelem += NAXIS1;
212 }
213 }
214
215
216 // Add the first binary table extension.
217 fits_create_tbl(fptr, BINARY_TBL, 1L, 3L, (char **)ttype1, (char **)tform1,
218 (char **)tunit1, NULL, &status);
219
220 // Write EXTNAME and EXTVER near the top, after TFIELDS.
221 fits_read_key_lng(fptr, "TFIELDS", &dummy, NULL, &status);
222 fits_insert_key_str(fptr, "EXTNAME", "WCS-TABLE",
223 "WCS Coordinate lookup table", &status);
224 fits_insert_key_lng(fptr, "EXTVER", 1L, "Table number 1", &status);
225
226 // Write the TDIM1 keyrecord after TFORM1.
227 fits_read_key_str(fptr, "TFORM1", keyrec, NULL, &status);
228 sprintf(keyrec, "(2,%ld,%ld)", K1, K2);
229 fits_insert_key_str(fptr, "TDIM1", keyrec, "Dimensions of 3-D array",
230 &status);
231
232 // Plate carr�e projection with a bit of noise for the sake of realism.
233 fp = array;
234 for (k2 = 0; k2 < K2; k2++) {
235 for (k1 = 0; k1 < K1; k1++) {
236 // Box-Muller transformation: uniform -> normal distribution.
237 x1 = lcprng();
238 x2 = lcprng();
239 if (x1 == 0.0) x1 = 1.0;
240 z = sqrt(-2.0 * log(x1));
241 x2 *= TWOPI;
242 z1 = z * cos(x2);
243 z2 = z * sin(x2);
244
245 *(fp++) = refval[0] + span[0] * (k1/(K1-1.0) - 0.5) + z1 * sigma[0];
246 *(fp++) = refval[1] + span[1] * (k2/(K2-1.0) - 0.5) + z2 * sigma[1];
247 }
248 }
249 fits_write_col_flt(fptr, 1, 1L, 1L, 2*K1*K2, array, &status);
250
251 fp = array;
252 for (k1 = 0; k1 < K1; k1++) {
253 *(fp++) = 4.0f * k1;
254 }
255 fits_write_col_flt(fptr, 2, 1L, 1L, K1, array, &status);
256
257 fp = array;
258 for (k2 = 0; k2 < K2; k2++) {
259 *(fp++) = 4.0f * k2;
260 }
261 fits_write_col_flt(fptr, 3, 1L, 1L, K2, array, &status);
262
263
264 // Add the second binary table extension.
265 if (fits_create_tbl(fptr, BINARY_TBL, 1L, 4L, (char **)ttype2,
266 (char **)tform2, (char **)tunit2, NULL, &status)) {
267 fits_report_error(stderr, status);
268 return 1;
269 }
270
271 // Write EXTNAME and EXTVER near the top, after TFIELDS.
272 fits_read_key_lng(fptr, "TFIELDS", &dummy, NULL, &status);
273 fits_insert_key_str(fptr, "EXTNAME", "WCS-TABLE",
274 "WCS Coordinate lookup table", &status);
275 fits_insert_key_lng(fptr, "EXTVER", 2L, "Table number 2", &status);
276
277 // Write the TDIM2 keyrecord after TFORM2.
278 fits_read_key_str(fptr, "TFORM2", keyrec, NULL, &status);
279 fits_insert_key_str(fptr, "TDIM2", "(1,8)", "Dimensions of 2-D array",
280 &status);
281
282 // Write the TDIM4 keyrecord after TFORM4.
283 fits_read_key_str(fptr, "TFORM4", keyrec, NULL, &status);
284 fits_insert_key_str(fptr, "TDIM4", "(1,8)", "Dimensions of 2-D array",
285 &status);
286
287
288 fits_write_col_dbl(fptr, 1, 1L, 1L, 8L, windex, &status);
289 fits_write_col_dbl(fptr, 2, 1L, 1L, 8L, wcoord, &status);
290 fits_write_col_dbl(fptr, 3, 1L, 1L, 8L, tindex, &status);
291 fits_write_col_dbl(fptr, 4, 1L, 1L, 8L, tcoord, &status);
292
293 fits_close_file(fptr, &status);
294
295 if (status) {
296 fits_report_error(stderr, status);
297 return 1;
298 }
299
300 return 0;
301 }
302
303 /*----------------------------------------------------------------------------
304 * A simple linear congruential pseudo-random number generator that produces
305 * the same results on all systems so that the test output can be compared.
306 * It produces a fixed sequence of uniformly distributed numbers in [0,1].
307 * Adapted from the example in Numerical Recipes in C.
308 *---------------------------------------------------------------------------*/
309
lcprng()310 double lcprng()
311 {
312 static unsigned long next = 137UL;
313
314 next = next * 1664525UL + 1013904223UL;
315 return (double)(next % 1073741824UL) / 1073741823.0;
316 }
317
318 //----------------------------------------------------------------------------
319
do_wcs_stuff(fitsfile * fptr,struct wcsprm * wcs)320 int do_wcs_stuff(fitsfile *fptr, struct wcsprm *wcs)
321
322 {
323 int i1, i2, i3, k, naxis1, naxis2, naxis3, stat[8], status;
324 double phi[8], pixcrd[8][4], imgcrd[8][4], theta[8], world[8][4],
325 x1, x2, x3;
326
327 // Initialize the wcsprm struct, also taking control of memory allocated by
328 // fits_read_wcstab().
329 if ((status = wcsset(wcs))) {
330 fprintf(stderr, "wcsset ERROR %d: %s.\n", status, wcs_errmsg[status]);
331 return 1;
332 }
333
334 // Print the struct.
335 if ((status = wcsprt(wcs))) return status;
336
337 // Compute coordinates in the corners.
338 fits_read_key(fptr, TINT, "NAXIS1", &naxis1, NULL, &status);
339 fits_read_key(fptr, TINT, "NAXIS2", &naxis2, NULL, &status);
340 fits_read_key(fptr, TINT, "NAXIS3", &naxis3, NULL, &status);
341
342 k = 0;
343 x3 = 1.0f;
344 for (i3 = 0; i3 < 2; i3++) {
345 x2 = 0.5f;
346
347 for (i2 = 0; i2 < 2; i2++) {
348 x1 = 0.5f;
349
350 for (i1 = 0; i1 < 2; i1++) {
351 pixcrd[k][0] = x1;
352 pixcrd[k][1] = x2;
353 pixcrd[k][2] = x3;
354 pixcrd[k][3] = 1.0f;
355
356 k++;
357 x1 = naxis1 + 0.5f;
358 }
359
360 x2 = naxis2 + 0.5f;
361 }
362
363 x3 = naxis3;
364 }
365
366 if ((status = wcsp2s(wcs, 8, 4, pixcrd[0], imgcrd[0], phi, theta, world[0],
367 stat))) {
368 fprintf(stderr, "\n\nwcsp2s ERROR %d: %s.\n", status,
369 wcs_errmsg[status]);
370
371 // Invalid pixel coordinates.
372 if (status == 8) status = 0;
373 }
374
375 if (status == 0) {
376 printf("\n\nCorner world coordinates:\n"
377 " Pixel World\n");
378 for (k = 0; k < 8; k++) {
379 printf(" (%5.1f,%6.1f,%4.1f,%4.1f) -> (%7.3f,%8.3f,%9g,%11.5f)",
380 pixcrd[k][0], pixcrd[k][1], pixcrd[k][2], pixcrd[k][3],
381 world[k][0], world[k][1], world[k][2], world[k][3]);
382 if (stat[k]) printf(" (BAD)");
383 printf("\n");
384 }
385 }
386
387 return 0;
388 }
389