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: tdis1.c,v 7.7 2021/07/12 06:36:49 mcalabre Exp $
23 *=============================================================================
24 *
25 * tdis1 tests the WCSLIB distortion functions for closure. Input comes from
26 * the FITS file specified as an argument, which must contain a 2D image, or
27 * else from TPV7.fits. The test is done via linp2x() and linx2p().
28 *
29 *---------------------------------------------------------------------------*/
30
31 #include <math.h>
32 #include <stdio.h>
33 #include <stdlib.h>
34 #include <string.h>
35 #include <time.h>
36
37 #include <wcserr.h>
38 #include <wcshdr.h>
39 #include <wcsprintf.h>
40 #include <lin.h>
41 #include <dis.h>
42
43 #define HEADER_SIZE 36000
44
45 int tpv2poly(struct wcsprm *wcs);
46
47 // Absolute and fractional tolerance. Distortions are typically used on
48 // large images, so the absolute tolerance in the corners may not be very
49 // high, simply due to floating point precision.
50 const double ATOL = 1.0e-9;
51 const double FTOL = 1.0e-10;
52
53
main(int argc,char * argv[])54 int main(int argc, char *argv[])
55
56 {
57 char *infile = "TPV7.fits";
58
59 char keyrec[81], header[288001], *disfn;
60 int dopoly, gotend, iblock, ikeyrec, inc, itest, j, k, n, naxis[2], naxis1,
61 naxis2, nClosure, nFail, nkeyrec, nsamp, nreject, nTest, nwcs, p1, p2,
62 status;
63 clock_t t0, tp2x, tx2p;
64 double absmax, dp1, dp2, *img, *img1, *img2, pix[8], pixblc[2], pixsamp[2],
65 pixtrc[2], px, *px0, *px1, pxi[8], rel, resid, relmax;
66 double *avgdis, *avgtot, *maxdis, *maxtot, *rmsdis, *rmstot, stats[9];
67 FILE *fptr;
68 struct linprm affine, *lin, *linpol, *lintpv;
69 struct wcsprm *wcs, wcspol;
70
71
72 wcserr_enable(1);
73 wcsprintf_set(stdout);
74
75 // Set line buffering in case stdout is redirected to a file, otherwise
76 // stdout and stderr messages will be jumbled (stderr is unbuffered).
77 setvbuf(stdout, NULL, _IOLBF, 0);
78
79 wcsprintf("Testing closure of WCSLIB distortion routines (tdis1.c)\n"
80 "-------------------------------------------------------\n");
81
82 // List status return messages.
83 wcsprintf("\nList of dis status return values:\n");
84 for (status = 1; status <= 5; status++) {
85 wcsprintf("%4d: %s.\n", status, dis_errmsg[status]);
86 }
87 wcsprintf("\n");
88
89 // Optional file name specified?
90 if (1 < argc) {
91 infile = argv[1];
92 }
93
94
95 // Read in the FITS header, excluding COMMENT and HISTORY keyrecords.
96 if ((fptr = fopen(infile, "r")) == 0) {
97 wcsprintf("ERROR opening %s\n", infile);
98 return 1;
99 }
100
101 memset(naxis, 0, 2*sizeof(int));
102
103 k = 0;
104 nkeyrec = 0;
105 gotend = 0;
106 for (iblock = 0; iblock < 100; iblock++) {
107 for (ikeyrec = 0; ikeyrec < 36; ikeyrec++) {
108 if (fgets(keyrec, 81, fptr) == 0) {
109 break;
110 }
111
112 if (strncmp(keyrec, " ", 8) == 0) continue;
113 if (strncmp(keyrec, "COMMENT ", 8) == 0) continue;
114 if (strncmp(keyrec, "HISTORY ", 8) == 0) continue;
115
116 if (strncmp(keyrec, "NAXIS", 5) == 0) {
117 if (keyrec[5] == ' ') {
118 sscanf(keyrec+10, "%d", &n);
119 if (n != 2) {
120 wcsprintf("ERROR, expecting a 2D image.\n");
121 return 1;
122 }
123 continue;
124 }
125
126 sscanf(keyrec+5, "%d = %d", &j, &n);
127 naxis[j-1] = n;
128 continue;
129 }
130
131 memcpy(header+k, keyrec, 80);
132 k += 80;
133 nkeyrec++;
134
135 if (strncmp(keyrec, "END ", 10) == 0) {
136 // An END keyrecord was read, but read the rest of the block.
137 gotend = 1;
138 }
139 }
140
141 if (gotend) break;
142 }
143 fclose(fptr);
144
145
146 // Parse the header.
147 if ((wcspih(header, nkeyrec, WCSHDR_none, 2, &nreject, &nwcs, &wcs))) {
148 wcsperr(wcs, 0x0);
149 return 1;
150 }
151
152 // Is it TPV?
153 dopoly = 0;
154 if (strcmp(wcs->ctype[0], "RA---TPV") == 0) {
155 // Copy it and translate to Polynomial for later use.
156 wcspol.flag = -1;
157 if (wcscopy(1, wcs, &wcspol)) {
158 wcsperr(wcs, 0x0);
159 return 1;
160 }
161
162 // Translate TPV to Polynomial.
163 tpv2poly(&wcspol);
164
165 wcspol.flag = -1;
166 if (wcsset(&wcspol)) {
167 wcsperr(&wcspol, 0x0);
168 return 1;
169 }
170
171 dopoly = 1;
172 }
173
174
175 // wcsset() translates the TPV "projection" into a sequent distortion.
176 if (wcsset(wcs)) {
177 wcsperr(wcs, 0x0);
178 return 1;
179 }
180
181 // Henceforth, we will work with linprm.
182 lin = &(wcs->lin);
183
184 // Get statistics on the distortion in the inner quarter of the image.
185 maxdis = stats;
186 maxtot = maxdis + 2;
187 avgdis = maxtot + 1;
188 avgtot = avgdis + 2;
189 rmsdis = avgtot + 1;
190 rmstot = rmsdis + 2;
191
192 pixblc[0] = 0.25 * naxis[0];
193 pixblc[1] = 0.25 * naxis[1];
194 pixtrc[0] = 0.75 * naxis[0];
195 pixtrc[1] = 0.75 * naxis[1];
196 pixsamp[0] = (pixtrc[0] - pixblc[0])/512.0;
197 pixsamp[1] = (pixtrc[1] - pixblc[1])/512.0;
198 if (pixsamp[0] < 1.0) pixsamp[0] = 1.0;
199 if (pixsamp[1] < 1.0) pixsamp[1] = 1.0;
200
201 if (linwarp(lin, pixblc, pixtrc, pixsamp, &nsamp,
202 maxdis, maxtot, avgdis, avgtot, rmsdis, rmstot)) {
203 linperr(lin, 0x0);
204 return 1;
205 }
206
207 for (k = 0; k < 9; k++) {
208 if (fabs(stats[k]) < 0.0005) stats[k] = 0.0;
209 }
210
211 wcsprintf("linwarp() statistics computed over %d sample points:\n"
212 " Max distortion, axis 1: %8.3f pixels\n"
213 " axis 2: %8.3f pixels\n"
214 " total: %8.3f pixels\n"
215 " Mean distortion, axis 1: %8.3f pixels\n"
216 " axis 2: %8.3f pixels\n"
217 " total: %8.3f pixels\n"
218 " RMS distortion, axis 1: %8.3f pixels\n"
219 " axis 2: %8.3f pixels\n"
220 " total: %8.3f pixels\n",
221 nsamp, maxdis[0], maxdis[1], *maxtot,
222 avgdis[0], avgdis[1], *avgtot,
223 rmsdis[0], rmsdis[1], *rmstot);
224
225 if (lin->disseq) {
226 // Exercise diswarp() as well.
227 wcsprintf("\n");
228
229 // Define a rectangle in intermediate pixel coordinates that just
230 // encompasses the inner quarter of the image. For this we need
231 // to switch off CDELTia scaling and all distortions.
232 affine.flag = -1;
233 if ((status = lincpy(1, lin, &affine))) {
234 linperr(lin, 0x0);
235 return 1;
236 }
237
238 affine.cdelt[0] = 1.0;
239 affine.cdelt[1] = 1.0;
240 if ((status = (lindis(1, &affine, 0x0) ||
241 lindis(2, &affine, 0x0) ||
242 linset(&affine)))) {
243 linperr(&affine, 0x0);
244 return 1;
245 }
246
247 pix[0] = pixblc[0];
248 pix[1] = pixblc[1];
249 pix[2] = pixtrc[0];
250 pix[3] = pixblc[1];
251 pix[4] = pixtrc[0];
252 pix[5] = pixtrc[1];
253 pix[6] = pixblc[0];
254 pix[7] = pixtrc[1];
255 if (linp2x(&affine, 4, 2, pix, pxi)) {
256 linperr(&affine, 0x0);
257 return 1;
258 }
259
260 linfree(&affine);
261
262 pixblc[0] = pxi[0];
263 pixblc[1] = pxi[1];
264 pixtrc[0] = pxi[0];
265 pixtrc[1] = pxi[1];
266 k = 2;
267 for (j = 1; j < 4; j++) {
268 if (pixblc[0] > pxi[k]) pixblc[0] = pxi[k];
269 if (pixtrc[0] < pxi[k]) pixtrc[0] = pxi[k];
270 k++;
271 if (pixblc[1] > pxi[k]) pixblc[1] = pxi[k];
272 if (pixtrc[1] < pxi[k]) pixtrc[1] = pxi[k];
273 k++;
274 }
275
276 pixsamp[0] = (pixtrc[0] - pixblc[0])/512.0;
277 pixsamp[1] = (pixtrc[1] - pixblc[1])/512.0;
278
279 if (diswarp(lin->disseq, pixblc, pixtrc, pixsamp, &nsamp,
280 maxdis, maxtot, avgdis, avgtot, rmsdis, rmstot)) {
281 wcserr_prt(lin->disseq->err, 0x0);
282 return 1;
283 }
284
285 for (k = 0; k < 9; k++) {
286 if (fabs(stats[k]) < 0.0005) stats[k] = 0.0;
287 }
288
289 wcsprintf("diswarp() statistics computed over %d sample points:\n"
290 " Max distortion, axis 1: %8.3f units\n"
291 " axis 2: %8.3f units\n"
292 " total: %8.3f units\n"
293 " Mean distortion, axis 1: %8.3f units\n"
294 " axis 2: %8.3f units\n"
295 " total: %8.3f units\n"
296 " RMS distortion, axis 1: %8.3f units\n"
297 " axis 2: %8.3f units\n"
298 " total: %8.3f units\n",
299 nsamp, maxdis[0], maxdis[1], *maxtot,
300 avgdis[0], avgdis[1], *avgtot,
301 rmsdis[0], rmsdis[1], *rmstot);
302 }
303
304
305 // The image size determines the test domain.
306 if ((naxis1 = naxis[0]) == 0) {
307 naxis1 = 2*wcs->crpix[0] + 1;
308 }
309 if ((naxis2 = naxis[1]) == 0) {
310 naxis2 = 2*wcs->crpix[1] + 1;
311 }
312
313 // Limit the number of tests.
314 inc = 1;
315 while ((naxis1/inc)*(naxis2/inc) > 800000) {
316 inc *= 2;
317 }
318
319 n = naxis1 / inc;
320 px0 = calloc(4*(2*n), sizeof(double));
321 px1 = px0 + 2*n ;
322 img = px1 + 2*n ;
323 img1 = img;
324 img2 = img + 2*n;
325
326 for (itest = 0; itest < 2; itest++) {
327 if (itest) {
328 if (!dopoly) break;
329
330 lin = &(wcspol.lin);
331 }
332
333 if (lin->dispre) {
334 disfn = lin->dispre->dtype[0];
335 } else if (lin->disseq) {
336 disfn = lin->disseq->dtype[0];
337 } else {
338 disfn = 0x0;
339 }
340
341 wcsprintf("\n");
342
343 // Now the closure test.
344 tp2x = 0;
345 tx2p = 0;
346 nTest = 0;
347 nFail = 0;
348 nClosure = 0;
349 absmax = 0.0;
350 relmax = 0.0;
351 for (p2 = 1; p2 <= naxis2; p2 += inc) {
352 k = 0;
353 for (p1 = 1; p1 <= naxis1; p1 += inc) {
354 px0[k++] = (double)p1;
355 px0[k++] = (double)p2;
356 }
357
358 t0 = clock();
359 if (linp2x(lin, n, 2, px0, img)) {
360 linperr(lin, 0x0);
361 nFail = 1;
362 break;
363 }
364 tp2x += clock() - t0;
365
366 t0 = clock();
367 if (linx2p(lin, n, 2, img, px1)) {
368 linperr(lin, 0x0);
369 nFail = 1;
370 break;
371 }
372 tx2p += clock() - t0;
373
374 // Check closure.
375 k = 0;
376 for (k = 0; k < 2*n ; k += 2) {
377 dp1 = fabs(px1[k] - px0[k]);
378 dp2 = fabs(px1[k+1] - px0[k+1]);
379
380 resid = (dp1 > dp2) ? dp1 : dp2;
381 if (resid > absmax) absmax = resid;
382
383 if (resid > ATOL) {
384 nClosure++;
385 wcsprintf("Absolute closure error:\n");
386 wcsprintf(" pix: %18.12f %18.12f\n", px0[k], px0[k+1]);
387 wcsprintf(" -> img: %18.12f %18.12f\n", img[k], img[k+1]);
388 wcsprintf(" -> pix: %18.12f %18.12f\n", px1[k], px1[k+1]);
389 wcsprintf("\n");
390 continue;
391 }
392
393 resid = 0.0;
394 if ((px = fabs(px0[k])) > 1.0) resid = dp1/px;
395 if ((px = fabs(px0[k+1])) > 1.0) {
396 if ((rel = dp2/px) > resid) resid = rel;
397 }
398 if (resid > relmax) relmax = resid;
399
400 if (resid > FTOL) {
401 nClosure++;
402 wcsprintf("Relative closure error:\n");
403 wcsprintf(" pix: %18.12f %18.12f\n", px0[k], px0[k+1]);
404 wcsprintf(" -> img: %18.12f %18.12f\n", img[k], img[k+1]);
405 wcsprintf(" -> pix: %18.12f %18.12f\n", px1[k], px1[k+1]);
406 wcsprintf("\n");
407 }
408 }
409
410 nTest += n;
411 }
412
413 if (nFail) {
414 wcsprintf("\nFAIL: The %s test failed to complete.\n", disfn);
415
416 } else {
417 wcsprintf("linp2x/linx2p with %s distortions:\n"
418 " Completed %d closure tests.\n"
419 " Maximum absolute closure residual = %.2e pixel.\n"
420 " Maximum relative closure residual = %.2e.\n", disfn,
421 nTest, absmax, relmax);
422 wcsprintf("\n");
423
424 wcsprintf(" linp2x time (ns): %6.0f\n linx2p time (ns): %6.0f\n\n",
425 1.0e9*((double)tp2x/CLOCKS_PER_SEC)/nTest,
426 1.0e9*((double)tx2p/CLOCKS_PER_SEC)/nTest);
427
428 if (nClosure) {
429 wcsprintf("FAIL: %d closure residuals exceed reporting tolerance.\n",
430 nClosure);
431
432 } else {
433 wcsprintf("PASS: All %s closure residuals are within reporting "
434 "tolerance.\n", disfn);
435 }
436 }
437 }
438
439
440 // Compare TPV with Polynomial over the test domain.
441 if (dopoly) {
442 wcsprintf("\n");
443
444 nTest = 0;
445 nFail = 0;
446 absmax = 0.0;
447 lintpv = &(wcs->lin);
448 linpol = &(wcspol.lin);
449 for (p2 = 1; p2 <= naxis2; p2 += inc) {
450 k = 0;
451 for (p1 = 1; p1 <= naxis1; p1 += inc) {
452 px0[k++] = (double)p1;
453 px0[k++] = (double)p2;
454 }
455
456 if (linp2x(lintpv, n, 2, px0, img1)) {
457 linperr(lintpv, 0x0);
458 break;
459 }
460
461 if (linp2x(linpol, n, 2, px0, img2)) {
462 linperr(linpol, 0x0);
463 break;
464 }
465
466 // Check agreement.
467 k = 0;
468 for (k = 0; k < 2*n ; k += 2) {
469 dp1 = fabs(img2[k] - img1[k]);
470 dp2 = fabs(img2[k+1] - img1[k+1]);
471
472 resid = (dp1 > dp2) ? dp1 : dp2;
473 if (resid > absmax) absmax = resid;
474
475 if (resid > ATOL) {
476 nFail++;
477 wcsprintf("TPV - Polynomial disagreement:\n");
478 wcsprintf(" pix: %18.12f %18.12f\n", px0[k], px0[k+1]);
479 wcsprintf(" -> TPV: %18.12f %18.12f\n", img1[k], img1[k+1]);
480 wcsprintf(" -> Pol: %18.12f %18.12f\n", img2[k], img2[k+1]);
481 wcsprintf("\n");
482 continue;
483 }
484 }
485
486 nTest += n;
487 }
488
489 wcsprintf("linp2x, TPV vs Polynomial distortions:\n"
490 " Completed %d comparisons.\n"
491 " Maximum absolute disagreement = %.2e units.\n", nTest, absmax);
492 wcsprintf("\n");
493
494 if (nFail) {
495 wcsprintf("FAIL: %d comparisons exceed reporting tolerance.\n", nFail);
496
497 } else {
498 wcsprintf("PASS: All TPV vs Polynomial comparisons are within "
499 "reporting tolerance.\n");
500 }
501 }
502
503
504 free(px0);
505 wcsvfree(&nwcs, &wcs);
506 wcsfree(&wcspol);
507
508 return nFail || nClosure;
509 }
510
511 /*----------------------------------------------------------------------------
512 * Translate a TPV "projection" to a general Polynomial distortion. Must be
513 * done before wcsset() which automatically translates PVi_ma keyvalues to a
514 * TPV distortion and elides them.
515 ----------------------------------------------------------------------------*/
516
tpv2poly(struct wcsprm * wcs)517 int tpv2poly(struct wcsprm *wcs)
518
519 {
520 const int tpvpow[40][3] = {
521 {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {2, 0, 0}, // 0 - 4
522 {1, 1, 0}, {0, 2, 0}, {3, 0, 0}, {2, 1, 0}, {1, 2, 0}, // 5 - 9
523 {0, 3, 0}, {0, 0, 3}, {4, 0, 0}, {3, 1, 0}, {2, 2, 0}, // 10 - 14
524 {1, 3, 0}, {0, 4, 0}, {5, 0, 0}, {4, 1, 0}, {3, 2, 0}, // 15 - 19
525 {2, 3, 0}, {1, 4, 0}, {0, 5, 0}, {0, 0, 5}, {6, 0, 0}, // 20 - 24
526 {5, 1, 0}, {4, 2, 0}, {3, 3, 0}, {2, 4, 0}, {1, 5, 0}, // 25 - 29
527 {0, 6, 0}, {7, 0, 0}, {6, 1, 0}, {5, 2, 0}, {4, 3, 0}, // 30 - 34
528 {3, 4, 0}, {2, 5, 0}, {1, 6, 0}, {0, 7, 0}, {0, 0, 7}}; // 35 - 39
529
530 char dqi[16], dqifield[32];
531 int i, k, m, naux, nterms;
532 struct disprm *dis;
533 struct dpkey *keyp;
534
535 // Attach a disprm struct to linprm. Also inits it.
536 if ((dis = calloc(1, sizeof(struct disprm))) == 0x0) {
537 return 1;
538 }
539
540 // Set values in the disprm struct.
541 dis->flag = -1;
542 disndp(24 + 3*wcs->npv);
543 lindis(2, &(wcs->lin), dis);
544
545
546 // The asterisk after the name prevents translation to TPD.
547 sprintf(dis->dtype[0], "Polynomial*");
548 sprintf(dis->dtype[1], "Polynomial*");
549
550 keyp = dis->dp;
551 for (i = 1; i <= 2; i++) {
552 // Glean essential information from the PVi_ma.
553 naux = 0;
554 nterms = 0;
555 for (k = 0; k < wcs->npv; k++) {
556 if (wcs->pv[k].i != i) continue;
557
558 m = wcs->pv[k].m;
559 if (tpvpow[m][2]) naux = 1;
560 nterms++;
561 }
562
563 // Distortion parameters for axis i.
564 sprintf(dqi, "DQ%d", i);
565 dpfill(keyp++, dqi, "NAXES", i, 0, 2, 0.0);
566 dpfill(keyp++, dqi, "AXIS.1", i, 0, (i==1)?1:2, 0.0);
567 dpfill(keyp++, dqi, "AXIS.2", i, 0, (i==1)?2:1, 0.0);
568 dpfill(keyp++, dqi, "DOCORR", i, 0, 0, 0.0);
569
570 if (naux) {
571 dpfill(keyp++, dqi, "NAUX", i, 0, 1, 0.0);
572 dpfill(keyp++, dqi, "AUX.1.COEFF.0", i, 1, 0, 0.0);
573 dpfill(keyp++, dqi, "AUX.1.POWER.0", i, 1, 0, 0.5);
574 dpfill(keyp++, dqi, "AUX.1.COEFF.1", i, 1, 0, 1.0);
575 dpfill(keyp++, dqi, "AUX.1.POWER.1", i, 1, 0, 2.0);
576 dpfill(keyp++, dqi, "AUX.1.COEFF.2", i, 1, 0, 1.0);
577 dpfill(keyp++, dqi, "AUX.1.POWER.2", i, 1, 0, 2.0);
578 }
579 dpfill(keyp++, dqi, "NTERMS", i, 0, nterms, 0.0);
580
581 nterms = 0;
582 for (k = 0; k < wcs->npv; k++) {
583 if (wcs->pv[k].i != i) continue;
584
585 sprintf(dqifield, "%s.TERM.%d", dqi, ++nterms);
586 dpfill(keyp++, dqifield, "COEFF", i, 1, 0, wcs->pv[k].value);
587
588 m = wcs->pv[k].m;
589 if (tpvpow[m][0]) {
590 dpfill(keyp++, dqifield, "VAR.1", i, 1, 0, (double)tpvpow[m][0]);
591 }
592
593 if (tpvpow[m][1]) {
594 dpfill(keyp++, dqifield, "VAR.2", i, 1, 0, (double)tpvpow[m][1]);
595 }
596
597 if (tpvpow[m][2]) {
598 dpfill(keyp++, dqifield, "AUX.1", i, 1, 0, (double)tpvpow[m][2]);
599 }
600 }
601 }
602
603 dis->ndp = keyp - dis->dp;
604
605
606 // TPV now becomes TAN.
607 strcpy(wcs->ctype[0]+5, "TAN");
608 strcpy(wcs->ctype[1]+5, "TAN");
609 wcs->npv = 0;
610
611 return 0;
612 }
613