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: wcsfix.c,v 7.7 2021/07/12 06:36:49 mcalabre Exp $
23 *===========================================================================*/
24
25 #include <math.h>
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <string.h>
29
30 #include "lin.h"
31 #include "sph.h"
32 #include "tab.h"
33 #include "wcs.h"
34 #include "wcserr.h"
35 #include "wcsfix.h"
36 #include "wcsmath.h"
37 #include "wcstrig.h"
38 #include "wcsunits.h"
39 #include "wcsutil.h"
40 #include "wtbarr.h"
41
42 extern const int WCSSET;
43
44 // Maximum number of coordinate axes that can be handled.
45 #define NMAX 16
46
47 // Map status return value to message.
48 const char *wcsfix_errmsg[] = {
49 "Success",
50 "Null wcsprm pointer passed",
51 "Memory allocation failed",
52 "Linear transformation matrix is singular",
53 "Inconsistent or unrecognized coordinate axis types",
54 "Invalid parameter value",
55 "Invalid coordinate transformation parameters",
56 "Ill-conditioned coordinate transformation parameters",
57 "All of the corner pixel coordinates are invalid",
58 "Could not determine reference pixel coordinate",
59 "Could not determine reference pixel value"};
60
61 // Map error returns for lower-level routines.
62 const int fix_linerr[] = {
63 FIXERR_SUCCESS, // 0: LINERR_SUCCESS
64 FIXERR_NULL_POINTER, // 1: LINERR_NULL_POINTER
65 FIXERR_MEMORY, // 2: LINERR_MEMORY
66 FIXERR_SINGULAR_MTX, // 3: LINERR_SINGULAR_MTX
67 FIXERR_BAD_PARAM, // 4: LINERR_DISTORT_INIT
68 FIXERR_NO_REF_PIX_COORD, // 5: LINERR_DISTORT
69 FIXERR_NO_REF_PIX_VAL // 6: LINERR_DEDISTORT
70 };
71
72 const int fix_wcserr[] = {
73 FIXERR_SUCCESS, // 0: WCSERR_SUCCESS
74 FIXERR_NULL_POINTER, // 1: WCSERR_NULL_POINTER
75 FIXERR_MEMORY, // 2: WCSERR_MEMORY
76 FIXERR_SINGULAR_MTX, // 3: WCSERR_SINGULAR_MTX
77 FIXERR_BAD_CTYPE, // 4: WCSERR_BAD_CTYPE
78 FIXERR_BAD_PARAM, // 5: WCSERR_BAD_PARAM
79 FIXERR_BAD_COORD_TRANS, // 6: WCSERR_BAD_COORD_TRANS
80 FIXERR_ILL_COORD_TRANS, // 7: WCSERR_ILL_COORD_TRANS
81 FIXERR_BAD_CORNER_PIX, // 8: WCSERR_BAD_PIX
82 FIXERR_NO_REF_PIX_VAL, // 9: WCSERR_BAD_WORLD
83 FIXERR_NO_REF_PIX_VAL // 10: WCSERR_BAD_WORLD_COORD
84 // ...others not used
85 };
86
87 // Convenience macro for invoking wcserr_set().
88 #define WCSFIX_ERRMSG(status) WCSERR_SET(status), wcsfix_errmsg[status]
89
90 //----------------------------------------------------------------------------
91
wcsfix(int ctrl,const int naxis[],struct wcsprm * wcs,int stat[])92 int wcsfix(int ctrl, const int naxis[], struct wcsprm *wcs, int stat[])
93
94 {
95 int status = 0;
96
97 if ((stat[CDFIX] = cdfix(wcs)) > 0) {
98 status = 1;
99 }
100
101 if ((stat[DATFIX] = datfix(wcs)) > 0) {
102 status = 1;
103 }
104
105 if ((stat[OBSFIX] = obsfix(0, wcs)) > 0) {
106 status = 1;
107 }
108
109 if ((stat[UNITFIX] = unitfix(ctrl, wcs)) > 0) {
110 status = 1;
111 }
112
113 if ((stat[SPCFIX] = spcfix(wcs)) > 0) {
114 status = 1;
115 }
116
117 if ((stat[CELFIX] = celfix(wcs)) > 0) {
118 status = 1;
119 }
120
121 if ((stat[CYLFIX] = cylfix(naxis, wcs)) > 0) {
122 status = 1;
123 }
124
125 return status;
126 }
127
128 //----------------------------------------------------------------------------
129
wcsfixi(int ctrl,const int naxis[],struct wcsprm * wcs,int stat[],struct wcserr info[])130 int wcsfixi(
131 int ctrl,
132 const int naxis[],
133 struct wcsprm *wcs,
134 int stat[],
135 struct wcserr info[])
136
137 {
138 int status = 0;
139
140 // Handling the status values returned from the sub-fixers is trickier than
141 // it might seem, especially considering that wcs->err may contain an error
142 // status on input which should be preserved if no translation errors occur.
143 // The simplest way seems to be to save a copy of wcs->err and clear it
144 // before each sub-fixer. The last real error to occur, excluding
145 // informative messages, is the one returned.
146
147 // To get informative messages from spcfix() it must precede celfix() and
148 // cylfix(). The latter call wcsset() which also translates AIPS-convention
149 // spectral axes.
150 struct wcserr err;
151 wcserr_copy(wcs->err, &err);
152
153 for (int ifix = CDFIX; ifix < NWCSFIX; ifix++) {
154 // Clear (delete) wcs->err.
155 wcserr_clear(&(wcs->err));
156
157 switch (ifix) {
158 case CDFIX:
159 stat[ifix] = cdfix(wcs);
160 break;
161 case DATFIX:
162 stat[ifix] = datfix(wcs);
163 break;
164 case OBSFIX:
165 stat[ifix] = obsfix(0, wcs);
166 break;
167 case UNITFIX:
168 stat[ifix] = unitfix(ctrl, wcs);
169 break;
170 case SPCFIX:
171 stat[ifix] = spcfix(wcs);
172 break;
173 case CELFIX:
174 stat[ifix] = celfix(wcs);
175 break;
176 case CYLFIX:
177 stat[ifix] = cylfix(naxis, wcs);
178 break;
179 default:
180 continue;
181 }
182
183 if (stat[ifix] == FIXERR_NO_CHANGE) {
184 // No change => no message.
185 wcserr_copy(0x0, info+ifix);
186
187 } else if (stat[ifix] == 0) {
188 // Successful translation, but there may be an informative message.
189 if (wcs->err && wcs->err->status < 0) {
190 wcserr_copy(wcs->err, info+ifix);
191 } else {
192 wcserr_copy(0x0, info+ifix);
193 }
194
195 } else {
196 // An informative message or error message.
197 wcserr_copy(wcs->err, info+ifix);
198
199 if ((status = (stat[ifix] > 0))) {
200 // It was an error, replace the previous one.
201 wcserr_copy(wcs->err, &err);
202 }
203 }
204 }
205
206 // Restore the last error to occur.
207 if (err.status) {
208 wcserr_copy(&err, wcs->err);
209 } else {
210 wcserr_clear(&(wcs->err));
211 }
212
213 return status;
214 }
215
216 //----------------------------------------------------------------------------
217
cdfix(struct wcsprm * wcs)218 int cdfix(struct wcsprm *wcs)
219
220 {
221 if (wcs == 0x0) return FIXERR_NULL_POINTER;
222
223 if ((wcs->altlin & 1) || !(wcs->altlin & 2)) {
224 // Either we have PCi_ja or there are no CDi_ja.
225 return FIXERR_NO_CHANGE;
226 }
227
228 int naxis = wcs->naxis;
229 int status = FIXERR_NO_CHANGE;
230 for (int i = 0; i < naxis; i++) {
231 // Row of zeros?
232 double *cd = wcs->cd + i*naxis;
233 for (int k = 0; k < naxis; k++, cd++) {
234 if (*cd != 0.0) goto next;
235 }
236
237 // Column of zeros?
238 cd = wcs->cd + i;
239 for (int k = 0; k < naxis; k++, cd += naxis) {
240 if (*cd != 0.0) goto next;
241 }
242
243 cd = wcs->cd + i * (naxis + 1);
244 *cd = 1.0;
245 status = FIXERR_SUCCESS;
246
247 next: ;
248 }
249
250 return status;
251 }
252
253 //----------------------------------------------------------------------------
254
parse_date(const char * buf,int * hour,int * minute,double * sec)255 static int parse_date(const char *buf, int *hour, int *minute, double *sec)
256
257 {
258 char ctmp[72];
259 if (sscanf(buf, "%2d:%2d:%s", hour, minute, ctmp) < 3 ||
260 wcsutil_str2double(ctmp, sec)) {
261 return 1;
262 }
263
264 return 0;
265 }
266
267
write_date(char * buf,int hour,int minute,double sec)268 static void write_date(char *buf, int hour, int minute, double sec)
269
270 {
271 char ctmp[32];
272 wcsutil_double2str(ctmp, "%04.1f", sec);
273 sprintf(buf, "T%.2d:%.2d:%s", hour, minute, ctmp);
274 }
275
276
newline(char ** cp)277 static char *newline(char **cp)
278
279 {
280 size_t k;
281 if ((k = strlen(*cp))) {
282 *cp += k;
283 strcat(*cp, ".\n");
284 *cp += 2;
285 }
286
287 return *cp;
288 }
289
290
datfix(struct wcsprm * wcs)291 int datfix(struct wcsprm *wcs)
292
293 {
294 static const char *function = "datfix";
295
296 // MJD of J2000.0 and B1900.0.
297 const double mjd2000 = 51544.5;
298 const double mjd1900 = 15019.81352;
299
300 // Days per Julian year and per tropical year.
301 const double djy = 365.25;
302 const double dty = 365.242198781;
303
304 int day, hour = 0, minute = 0, month, year;
305 double sec = 0.0;
306
307 if (wcs == 0x0) return FIXERR_NULL_POINTER;
308 struct wcserr **err = &(wcs->err);
309
310 char infomsg[512];
311 char *cp = infomsg;
312 *cp = '\0';
313
314 int status = FIXERR_NO_CHANGE;
315
316 for (int i = 0; i < 5; i++) {
317 // MJDREF is split into integer and fractional parts, wheres MJDOBS and
318 // the rest are a single value.
319 const char *dateid;
320 char *date;
321 double *wcsmjd;
322 if (i == 0) {
323 // Note, DATEREF and MJDREF, not DATE-REF and MJD-REF (sigh).
324 dateid = "REF";
325 date = wcs->dateref;
326 wcsmjd = wcs->mjdref;
327 } else if (i == 1) {
328 dateid = "-OBS";
329 date = wcs->dateobs;
330 wcsmjd = &(wcs->mjdobs);
331 } else if (i == 2) {
332 dateid = "-BEG";
333 date = wcs->datebeg;
334 wcsmjd = &(wcs->mjdbeg);
335 } else if (i == 3) {
336 dateid = "-AVG";
337 date = wcs->dateavg;
338 wcsmjd = &(wcs->mjdavg);
339 } else if (i == 4) {
340 dateid = "-END";
341 date = wcs->dateend;
342 wcsmjd = &(wcs->mjdend);
343 }
344
345 char orig_date[72];
346 strncpy(orig_date, date, 72);
347
348 if (date[0] == '\0') {
349 // Fill in DATE from MJD if possible.
350
351 if (i == 1 && undefined(*wcsmjd)) {
352 // See if we have jepoch or bepoch.
353 if (!undefined(wcs->jepoch)) {
354 *wcsmjd = mjd2000 + (wcs->jepoch - 2000.0)*djy;
355 sprintf(newline(&cp), "Set MJD-OBS to %.6f from JEPOCH", *wcsmjd);
356 if (status == FIXERR_NO_CHANGE) status = FIXERR_SUCCESS;
357
358 } else if (!undefined(wcs->bepoch)) {
359 *wcsmjd = mjd1900 + (wcs->bepoch - 1900.0)*dty;
360 sprintf(newline(&cp), "Set MJD-OBS to %.6f from BEPOCH", *wcsmjd);
361 if (status == FIXERR_NO_CHANGE) status = FIXERR_SUCCESS;
362 }
363 }
364
365 if (undefined(*wcsmjd)) {
366 // No date information was provided.
367
368 } else {
369 // Calendar date from MJD, with allowance for MJD < 0.
370 double mjd[2], t;
371 if (i == 0) {
372 // MJDREF is already split into integer and fractional parts.
373 mjd[0] = wcsmjd[0];
374 mjd[1] = wcsmjd[1];
375 if (1.0 < mjd[1]) {
376 // Ensure the fractional part lies between 0 and +1.
377 t = floor(mjd[1]);
378 mjd[0] += t;
379 mjd[1] -= t;
380 }
381 } else {
382 // Split it into integer and fractional parts.
383 mjd[0] = floor(*wcsmjd);
384 mjd[1] = *wcsmjd - mjd[0];
385 }
386
387 int jd = 2400001 + (int)mjd[0];
388
389 int n4 = 4*(jd + ((2*((4*jd - 17918)/146097)*3)/4 + 1)/2 - 37);
390 int dd = 10*(((n4-237)%1461)/4) + 5;
391
392 year = n4/1461 - 4712;
393 month = (2 + dd/306)%12 + 1;
394 day = (dd%306)/10 + 1;
395 sprintf(date, "%.4d-%.2d-%.2d", year, month, day);
396
397 // Write time part only if non-zero.
398 if (0.0 < (t = mjd[1])) {
399 t *= 24.0;
400 hour = (int)t;
401 t = 60.0 * (t - hour);
402 minute = (int)t;
403 sec = 60.0 * (t - minute);
404
405 // Round to 1ms.
406 dd = 60000*(60*hour + minute) + (int)(1000*(sec+0.0005));
407 hour = dd / 3600000;
408 dd -= 3600000 * hour;
409 minute = dd / 60000;
410 int msec = dd - 60000 * minute;
411 sprintf(date+10, "T%.2d:%.2d:%.2d", hour, minute, msec/1000);
412
413 // Write fractions of a second only if non-zero.
414 if (msec%1000) {
415 sprintf(date+19, ".%.3d", msec%1000);
416 }
417 }
418 }
419
420 } else {
421 if (strlen(date) < 8) {
422 // Can't be a valid date.
423 status = FIXERR_BAD_PARAM;
424 sprintf(newline(&cp), "Invalid DATE%s format '%s' is too short",
425 dateid, date);
426 continue;
427 }
428
429 // Identify the date format.
430 if (date[4] == '-' && date[7] == '-') {
431 // Standard year-2000 form: CCYY-MM-DD[Thh:mm:ss[.sss...]]
432 if (sscanf(date, "%4d-%2d-%2d", &year, &month, &day) < 3) {
433 status = FIXERR_BAD_PARAM;
434 sprintf(newline(&cp), "Invalid DATE%s format '%s'", dateid, date);
435 continue;
436 }
437
438 if (date[10] == 'T') {
439 if (parse_date(date+11, &hour, &minute, &sec)) {
440 status = FIXERR_BAD_PARAM;
441 sprintf(newline(&cp), "Invalid time in DATE%s '%s'", dateid,
442 date+11);
443 continue;
444 }
445 } else if (date[10] == ' ') {
446 hour = 0;
447 minute = 0;
448 sec = 0.0;
449 if (parse_date(date+11, &hour, &minute, &sec)) {
450 write_date(date+10, hour, minute, sec);
451 } else {
452 date[10] = 'T';
453 }
454 }
455
456 } else if (date[4] == '/' && date[7] == '/') {
457 // Also allow CCYY/MM/DD[Thh:mm:ss[.sss...]]
458 if (sscanf(date, "%4d/%2d/%2d", &year, &month, &day) < 3) {
459 status = FIXERR_BAD_PARAM;
460 sprintf(newline(&cp), "Invalid DATE%s format '%s'", dateid, date);
461 continue;
462 }
463
464 if (date[10] == 'T') {
465 if (parse_date(date+11, &hour, &minute, &sec)) {
466 status = FIXERR_BAD_PARAM;
467 sprintf(newline(&cp), "Invalid time in DATE%s '%s'", dateid,
468 date+11);
469 continue;
470 }
471 } else if (date[10] == ' ') {
472 hour = 0;
473 minute = 0;
474 sec = 0.0;
475 if (parse_date(date+11, &hour, &minute, &sec)) {
476 write_date(date+10, hour, minute, sec);
477 } else {
478 date[10] = 'T';
479 }
480 }
481
482 // Looks ok, fix it up.
483 date[4] = '-';
484 date[7] = '-';
485
486 } else {
487 if (i == 1 && date[2] == '/' && date[5] == '/') {
488 // Old format DATE-OBS date: DD/MM/YY, also allowing DD/MM/CCYY.
489 if (sscanf(date, "%2d/%2d/%4d", &day, &month, &year) < 3) {
490 status = FIXERR_BAD_PARAM;
491 sprintf(newline(&cp), "Invalid DATE%s format '%s'", dateid,
492 date);
493 continue;
494 }
495
496 } else if (i == 1 && date[2] == '-' && date[5] == '-') {
497 // Also recognize DD-MM-YY and DD-MM-CCYY
498 if (sscanf(date, "%2d-%2d-%4d", &day, &month, &year) < 3) {
499 status = FIXERR_BAD_PARAM;
500 sprintf(newline(&cp), "Invalid DATE%s format '%s'", dateid,
501 date);
502 continue;
503 }
504
505 } else {
506 // Not a valid date format.
507 status = FIXERR_BAD_PARAM;
508 sprintf(newline(&cp), "Invalid DATE%s format '%s'", dateid, date);
509 continue;
510 }
511
512 if (year < 100) year += 1900;
513
514 // Doesn't have a time.
515 sprintf(date, "%.4d-%.2d-%.2d", year, month, day);
516 }
517
518 // Compute MJD.
519 double mjd[2];
520 mjd[0] = (double)((1461*(year - (12-month)/10 + 4712))/4
521 + (306*((month+9)%12) + 5)/10
522 - (3*((year - (12-month)/10 + 4900)/100))/4
523 + day - 2399904);
524 mjd[1] = (hour + (minute + sec/60.0)/60.0)/24.0;
525 double mjdsum = mjd[0] + mjd[1];
526
527 if (undefined(*wcsmjd)) {
528 if (i == 0) {
529 wcsmjd[0] = mjd[0];
530 wcsmjd[1] = mjd[1];
531 } else {
532 *wcsmjd = mjdsum;
533 }
534 sprintf(newline(&cp), "Set MJD%s to %.6f from DATE%s", dateid,
535 mjdsum, dateid);
536
537 if (status == FIXERR_NO_CHANGE) status = FIXERR_SUCCESS;
538
539 } else {
540 // Check for consistency.
541 double mjdtmp;
542 if (i == 0) {
543 mjdtmp = wcsmjd[0] + wcsmjd[1];
544 } else {
545 mjdtmp = *wcsmjd;
546 }
547
548 if (0.001 < fabs(mjdsum - mjdtmp)) {
549 status = FIXERR_BAD_PARAM;
550 sprintf(newline(&cp),
551 "Invalid parameter values: MJD%s and DATE%s are inconsistent",
552 dateid, dateid);
553 }
554 }
555
556 if (i == 1) {
557 if (!undefined(wcs->jepoch)) {
558 // Check consistency of JEPOCH.
559 double jepoch = 2000.0 + (*wcsmjd - mjd2000) / djy;
560
561 if (0.000002 < fabs(jepoch - wcs->jepoch)) {
562 // Informational only, no error.
563 sprintf(newline(&cp), "JEPOCH is inconsistent with DATE-OBS");
564 }
565 }
566
567 if (!undefined(wcs->bepoch)) {
568 // Check consistency of BEPOCH.
569 double bepoch = 1900.0 + (*wcsmjd - mjd1900) / dty;
570
571 if (0.000002 < fabs(bepoch - wcs->bepoch)) {
572 // Informational only, no error.
573 sprintf(newline(&cp), "BEPOCH is inconsistent with DATE-OBS");
574 }
575 }
576 }
577 }
578
579 if (strncmp(orig_date, date, 72)) {
580 if (orig_date[0] == '\0') {
581 sprintf(newline(&cp), "Set DATE%s to '%s' from MJD%s", dateid, date,
582 dateid);
583 } else {
584 sprintf(newline(&cp), "Changed DATE%s from '%s' to '%s'", dateid,
585 orig_date, date);
586 }
587
588 if (status == FIXERR_NO_CHANGE) status = FIXERR_SUCCESS;
589 }
590 }
591
592 if (*infomsg) {
593 wcserr_set(WCSERR_SET(FIXERR_DATE_FIX), infomsg);
594 }
595
596 return status;
597 }
598
599 //----------------------------------------------------------------------------
600
obsfix(int ctrl,struct wcsprm * wcs)601 int obsfix(int ctrl, struct wcsprm *wcs)
602
603 {
604 static const char *function = "obsfix";
605
606 // IAU(1976) ellipsoid (as prescribed by WCS Paper VII).
607 const double a = 6378140.0;
608 const double f = 1.0 / 298.2577;
609 const double e2 = (2.0 - f)*f;
610
611 if (wcs == 0x0) return FIXERR_NULL_POINTER;
612 struct wcserr **err = &(wcs->err);
613
614 // Set masks for checking partially-defined coordinate triplets.
615 int havexyz = 7;
616 havexyz -= 1*undefined(wcs->obsgeo[0]);
617 havexyz -= 2*undefined(wcs->obsgeo[1]);
618 havexyz -= 4*undefined(wcs->obsgeo[2]);
619 int havelbh = 7;
620 havelbh -= 1*undefined(wcs->obsgeo[3]);
621 havelbh -= 2*undefined(wcs->obsgeo[4]);
622 havelbh -= 4*undefined(wcs->obsgeo[5]);
623
624 if (ctrl == 2) {
625 // Make no changes.
626 if (0 < havexyz && havexyz < 7) {
627 return wcserr_set(WCSERR_SET(FIXERR_BAD_PARAM),
628 "Partially undefined Cartesian coordinate triplet");
629 }
630
631 if (0 < havelbh && havelbh < 7) {
632 return wcserr_set(WCSERR_SET(FIXERR_BAD_PARAM),
633 "Partially undefined Geodetic coordinate triplet");
634 }
635
636 if (havexyz == 0 || havelbh == 0) {
637 return FIXERR_NO_CHANGE;
638 }
639 }
640
641 if (havexyz == 0 && havelbh == 0) {
642 return FIXERR_NO_CHANGE;
643 }
644
645
646 char infomsg[256];
647 infomsg[0] = '\0';
648
649 int status = FIXERR_NO_CHANGE;
650
651 size_t k;
652 double x, y, z;
653 if (havelbh == 7) {
654 // Compute (x,y,z) from (lng,lat,hgt).
655 double coslat, coslng, sinlat, sinlng;
656 sincosd(wcs->obsgeo[3], &sinlng, &coslng);
657 sincosd(wcs->obsgeo[4], &sinlat, &coslat);
658 double n = a / sqrt(1.0 - e2*sinlat*sinlat);
659 double rho = n + wcs->obsgeo[5];
660
661 x = rho*coslng*coslat;
662 y = rho*sinlng*coslat;
663 z = (rho - n*e2)*sinlat;
664
665 if (havexyz < 7) {
666 // One or more of the Cartesian elements was undefined.
667 status = FIXERR_SUCCESS;
668 char *cp = infomsg;
669
670 if (ctrl == 1 || !(havexyz & 1)) {
671 wcs->obsgeo[0] = x;
672 sprintf(cp, "%s OBSGEO-X to %12.3f from OBSGEO-[LBH]",
673 (havexyz & 1) ? "Reset" : "Set", x);
674 }
675
676 if (ctrl == 1 || !(havexyz & 2)) {
677 wcs->obsgeo[1] = y;
678
679 if ((k = strlen(cp))) {
680 strcat(cp+k, ".\n");
681 cp += k + 2;
682 }
683
684 sprintf(cp, "%s OBSGEO-Y to %12.3f from OBSGEO-[LBH]",
685 (havexyz & 2) ? "Reset" : "Set", y);
686 }
687
688 if (ctrl == 1 || !(havexyz & 4)) {
689 wcs->obsgeo[2] = z;
690 if ((k = strlen(cp))) {
691 strcat(cp+k, ".\n");
692 cp += k + 2;
693 }
694
695 sprintf(cp, "%s OBSGEO-Z to %12.3f from OBSGEO-[LBH]",
696 (havexyz & 4) ? "Reset" : "Set", z);
697 }
698
699 wcserr_set(WCSERR_SET(FIXERR_OBSGEO_FIX), infomsg);
700
701 if (havexyz == 0) {
702 // Skip the consistency check.
703 return status;
704 }
705 }
706
707 } else if (havexyz == 7) {
708 // Compute (lng,lat,hgt) from (x,y,z).
709 x = wcs->obsgeo[0];
710 y = wcs->obsgeo[1];
711 z = wcs->obsgeo[2];
712 double r2 = x*x + y*y;
713
714 // Iterate over the value of zeta.
715 double coslat, coslng, sinlat, sinlng;
716 double n, rho, zeta = z;
717 for (int i = 0; i < 4; i++) {
718 rho = sqrt(r2 + zeta*zeta);
719 sinlat = zeta / rho;
720 n = a / sqrt(1.0 - e2*sinlat*sinlat);
721
722 zeta = z / (1.0 - n*e2/rho);
723 }
724
725 double lng = atan2d(y, x);
726 double lat = asind(sinlat);
727 double hgt = rho - n;
728
729 if (havelbh < 7) {
730 // One or more of the Geodetic elements was undefined.
731 status = FIXERR_SUCCESS;
732 char *cp = infomsg;
733
734 if (ctrl == 1 || !(havelbh & 1)) {
735 wcs->obsgeo[3] = lng;
736 sprintf(cp, "%s OBSGEO-L to %12.6f from OBSGEO-[XYZ]",
737 (havelbh & 1) ? "Reset" : "Set", lng);
738 }
739
740 if (ctrl == 1 || !(havelbh & 2)) {
741 wcs->obsgeo[4] = lat;
742 if ((k = strlen(cp))) {
743 strcat(cp+k, ".\n");
744 cp += k + 2;
745 }
746
747 sprintf(cp, "%s OBSGEO-B to %12.6f from OBSGEO-[XYZ]",
748 (havelbh & 2) ? "Reset" : "Set", lat);
749 }
750
751 if (ctrl == 1 || !(havelbh & 4)) {
752 wcs->obsgeo[5] = hgt;
753 if ((k = strlen(cp))) {
754 strcat(cp+k, ".\n");
755 cp += k + 2;
756 }
757
758 sprintf(cp, "%s OBSGEO-H to %12.3f from OBSGEO-[XYZ]",
759 (havelbh & 4) ? "Reset" : "Set", hgt);
760 }
761
762 wcserr_set(WCSERR_SET(FIXERR_OBSGEO_FIX), infomsg);
763
764 if (havelbh == 0) {
765 // Skip the consistency check.
766 return status;
767 }
768 }
769
770 // Compute (x,y,z) from (lng,lat,hgt) for consistency checking.
771 sincosd(wcs->obsgeo[3], &sinlng, &coslng);
772 sincosd(wcs->obsgeo[4], &sinlat, &coslat);
773 n = a / sqrt(1.0 - e2*sinlat*sinlat);
774 rho = n + wcs->obsgeo[5];
775
776 x = rho*coslng*coslat;
777 y = rho*sinlng*coslat;
778 z = (rho - n*e2)*sinlat;
779
780 } else {
781 return wcserr_set(WCSERR_SET(FIXERR_BAD_PARAM),
782 "Observatory coordinates incomplete");
783 }
784
785
786 // Check consistency.
787 double d, r2 = 0.0;
788 d = wcs->obsgeo[0] - x;
789 r2 += d*d;
790 d = wcs->obsgeo[1] - y;
791 r2 += d*d;
792 d = wcs->obsgeo[2] - z;
793 r2 += d*d;
794
795 if (1.0 < r2) {
796 d = sqrt(r2);
797 return wcserr_set(WCSERR_SET(FIXERR_BAD_PARAM),
798 "Observatory coordinates inconsistent by %.1f metres", d);
799 }
800
801 return status;
802 }
803
804
805 //----------------------------------------------------------------------------
806
unitfix(int ctrl,struct wcsprm * wcs)807 int unitfix(int ctrl, struct wcsprm *wcs)
808
809 {
810 const char *function = "unitfix";
811
812 if (wcs == 0x0) return FIXERR_NULL_POINTER;
813 struct wcserr **err = &(wcs->err);
814
815 int status = FIXERR_NO_CHANGE;
816
817 char msg[512];
818 strncpy(msg, "Changed units:", 512);
819
820 for (int i = 0; i < wcs->naxis; i++) {
821 char orig_unit[72];
822 strncpy(orig_unit, wcs->cunit[i], 71);
823 int result = wcsutrne(ctrl, wcs->cunit[i], &(wcs->err));
824 if (result == 0 || result == 12) {
825 size_t msglen = strlen(msg);
826 if (msglen < 511) {
827 wcsutil_null_fill(72, orig_unit);
828 char msgtmp[192];
829 sprintf(msgtmp, "\n '%s' -> '%s',", orig_unit, wcs->cunit[i]);
830 strncpy(msg+msglen, msgtmp, 511-msglen);
831 status = FIXERR_UNITS_ALIAS;
832 }
833 }
834 }
835
836 if (status == FIXERR_UNITS_ALIAS) {
837 // Chop off the trailing ", ".
838 size_t msglen = strlen(msg) - 2;
839 msg[msglen] = '\0';
840 wcserr_set(WCSERR_SET(FIXERR_UNITS_ALIAS), msg);
841
842 status = FIXERR_SUCCESS;
843 }
844
845 return status;
846 }
847
848 //----------------------------------------------------------------------------
849
spcfix(struct wcsprm * wcs)850 int spcfix(struct wcsprm *wcs)
851
852 {
853 static const char *function = "spcfix";
854
855 if (wcs == 0x0) return FIXERR_NULL_POINTER;
856 struct wcserr **err = &(wcs->err);
857
858 for (int i = 0; i < wcs->naxis; i++) {
859 // Translate an AIPS-convention spectral type if present.
860 char ctype[9], specsys[9];
861 int status = spcaips(wcs->ctype[i], wcs->velref, ctype, specsys);
862 if (status == FIXERR_SUCCESS) {
863 // An AIPS type was found but it may match what we already have.
864 status = FIXERR_NO_CHANGE;
865
866 // Was specsys translated?
867 if (wcs->specsys[0] == '\0' && *specsys) {
868 strncpy(wcs->specsys, specsys, 9);
869 wcserr_set(WCSERR_SET(FIXERR_SPC_UPDATE),
870 "Changed SPECSYS to '%s'", specsys);
871 status = FIXERR_SUCCESS;
872 }
873
874 // Was ctype translated? Have to null-fill for comparing them.
875 wcsutil_null_fill(9, wcs->ctype[i]);
876 if (strncmp(wcs->ctype[i], ctype, 9)) {
877 // ctype was translated...
878 if (status == FIXERR_SUCCESS) {
879 // ...and specsys was also.
880 wcserr_set(WCSERR_SET(FIXERR_SPC_UPDATE),
881 "Changed CTYPE%d from '%s' to '%s', and SPECSYS to '%s' "
882 "(VELREF=%d)", i+1, wcs->ctype[i], ctype, wcs->specsys,
883 wcs->velref);
884 } else {
885 wcserr_set(WCSERR_SET(FIXERR_SPC_UPDATE),
886 "Changed CTYPE%d from '%s' to '%s' (VELREF=%d)", i+1,
887 wcs->ctype[i], ctype, wcs->velref);
888 status = FIXERR_SUCCESS;
889 }
890
891 strncpy(wcs->ctype[i], ctype, 9);
892 }
893
894 // Tidy up.
895 if (status == FIXERR_SUCCESS) {
896 wcsutil_null_fill(72, wcs->ctype[i]);
897 wcsutil_null_fill(72, wcs->specsys);
898 }
899
900 // No need to check for others, wcsset() will fail if so.
901 return status;
902
903 } else if (status == SPCERR_BAD_SPEC_PARAMS) {
904 // An AIPS spectral type was found but with invalid velref.
905 return wcserr_set(WCSERR_SET(FIXERR_BAD_PARAM),
906 "Invalid parameter value: velref = %d", wcs->velref);
907 }
908 }
909
910 return FIXERR_NO_CHANGE;
911 }
912
913 //----------------------------------------------------------------------------
914
celfix(struct wcsprm * wcs)915 int celfix(struct wcsprm *wcs)
916
917 {
918 static const char *function = "celfix";
919
920 if (wcs == 0x0) return FIXERR_NULL_POINTER;
921 struct wcserr **err = &(wcs->err);
922
923 // Initialize if required.
924 int status;
925 if (wcs->flag != WCSSET) {
926 if ((status = wcsset(wcs))) return fix_wcserr[status];
927 }
928
929 // Was an NCP or GLS projection code translated?
930 if (wcs->lat >= 0) {
931 // Check ctype.
932 if (strcmp(wcs->ctype[wcs->lat]+5, "NCP") == 0) {
933 strcpy(wcs->ctype[wcs->lng]+5, "SIN");
934 strcpy(wcs->ctype[wcs->lat]+5, "SIN");
935
936 if (wcs->npvmax < wcs->npv + 2) {
937 // Allocate space for two more PVi_ma keyvalues.
938 if (wcs->m_flag == WCSSET && wcs->pv == wcs->m_pv) {
939 if (!(wcs->pv = calloc(wcs->npv+2, sizeof(struct pvcard)))) {
940 wcs->pv = wcs->m_pv;
941 return wcserr_set(WCSFIX_ERRMSG(FIXERR_MEMORY));
942 }
943
944 wcs->npvmax = wcs->npv + 2;
945 wcs->m_flag = WCSSET;
946
947 for (int k = 0; k < wcs->npv; k++) {
948 wcs->pv[k] = wcs->m_pv[k];
949 }
950
951 if (wcs->m_pv) free(wcs->m_pv);
952 wcs->m_pv = wcs->pv;
953
954 } else {
955 return wcserr_set(WCSFIX_ERRMSG(FIXERR_MEMORY));
956 }
957 }
958
959 struct celprm *wcscel = &(wcs->cel);
960 struct prjprm *wcsprj = &(wcscel->prj);
961 wcs->pv[wcs->npv].i = wcs->lat + 1;
962 wcs->pv[wcs->npv].m = 1;
963 wcs->pv[wcs->npv].value = wcsprj->pv[1];
964 (wcs->npv)++;
965
966 wcs->pv[wcs->npv].i = wcs->lat + 1;
967 wcs->pv[wcs->npv].m = 2;
968 wcs->pv[wcs->npv].value = wcsprj->pv[2];
969 (wcs->npv)++;
970
971 return FIXERR_SUCCESS;
972
973 } else if (strcmp(wcs->ctype[wcs->lat]+5, "GLS") == 0) {
974 strcpy(wcs->ctype[wcs->lng]+5, "SFL");
975 strcpy(wcs->ctype[wcs->lat]+5, "SFL");
976
977 if (wcs->crval[wcs->lng] != 0.0 || wcs->crval[wcs->lat] != 0.0) {
978 // In the AIPS convention, setting the reference longitude and
979 // latitude for GLS does not create an oblique graticule. A non-zero
980 // reference longitude introduces an offset in longitude in the normal
981 // way, whereas a non-zero reference latitude simply translates the
982 // reference point (i.e. the map as a whole) to that latitude. This
983 // might be effected by adjusting CRPIXja but that is complicated by
984 // the linear transformation and instead is accomplished here by
985 // setting theta_0.
986 if (wcs->npvmax < wcs->npv + 3) {
987 // Allocate space for three more PVi_ma keyvalues.
988 if (wcs->m_flag == WCSSET && wcs->pv == wcs->m_pv) {
989 if (!(wcs->pv = calloc(wcs->npv+3, sizeof(struct pvcard)))) {
990 wcs->pv = wcs->m_pv;
991 return wcserr_set(WCSFIX_ERRMSG(FIXERR_MEMORY));
992 }
993
994 wcs->npvmax = wcs->npv + 3;
995 wcs->m_flag = WCSSET;
996
997 for (int k = 0; k < wcs->npv; k++) {
998 wcs->pv[k] = wcs->m_pv[k];
999 }
1000
1001 if (wcs->m_pv) free(wcs->m_pv);
1002 wcs->m_pv = wcs->pv;
1003
1004 } else {
1005 return wcserr_set(WCSFIX_ERRMSG(FIXERR_MEMORY));
1006 }
1007 }
1008
1009 wcs->pv[wcs->npv].i = wcs->lng + 1;
1010 wcs->pv[wcs->npv].m = 0;
1011 wcs->pv[wcs->npv].value = 1.0;
1012 (wcs->npv)++;
1013
1014 // Note that the reference longitude is still zero.
1015 wcs->pv[wcs->npv].i = wcs->lng + 1;
1016 wcs->pv[wcs->npv].m = 1;
1017 wcs->pv[wcs->npv].value = 0.0;
1018 (wcs->npv)++;
1019
1020 wcs->pv[wcs->npv].i = wcs->lng + 1;
1021 wcs->pv[wcs->npv].m = 2;
1022 wcs->pv[wcs->npv].value = wcs->crval[wcs->lat];
1023 (wcs->npv)++;
1024 }
1025
1026 return FIXERR_SUCCESS;
1027 }
1028 }
1029
1030 return FIXERR_NO_CHANGE;
1031 }
1032
1033 //----------------------------------------------------------------------------
1034
cylfix(const int naxis[],struct wcsprm * wcs)1035 int cylfix(const int naxis[], struct wcsprm *wcs)
1036
1037 {
1038 static const char *function = "cylfix";
1039
1040 if (naxis == 0x0) return FIXERR_NO_CHANGE;
1041 if (wcs == 0x0) return FIXERR_NULL_POINTER;
1042 struct wcserr **err = &(wcs->err);
1043
1044 // Initialize if required.
1045 int status;
1046 if (wcs->flag != WCSSET) {
1047 if ((status = wcsset(wcs))) return fix_wcserr[status];
1048 }
1049
1050 // Check that we have a cylindrical projection.
1051 if (wcs->cel.prj.category != CYLINDRICAL) return FIXERR_NO_CHANGE;
1052 if (wcs->naxis < 2) return FIXERR_NO_CHANGE;
1053
1054
1055 // Compute the native longitude in each corner of the image.
1056 unsigned short ncnr = 1 << wcs->naxis;
1057
1058 unsigned short indx[NMAX];
1059 for (int k = 0; k < NMAX; k++) {
1060 indx[k] = 1 << k;
1061 }
1062
1063 int stat[4];
1064 double img[4][NMAX], phi[4], pix[4][NMAX], theta[4], world[4][NMAX];
1065
1066 double phimin = 1.0e99;
1067 double phimax = -1.0e99;
1068 for (unsigned short icnr = 0; icnr < ncnr;) {
1069 // Do four corners at a time.
1070 for (int j = 0; j < 4; j++, icnr++) {
1071 double *pixj = pix[j];
1072
1073 for (int k = 0; k < wcs->naxis; k++) {
1074 if (icnr & indx[k]) {
1075 *(pixj++) = naxis[k] + 0.5;
1076 } else {
1077 *(pixj++) = 0.5;
1078 }
1079 }
1080 }
1081
1082 if (!(status = wcsp2s(wcs, 4, NMAX, pix[0], img[0], phi, theta, world[0],
1083 stat))) {
1084 for (int j = 0; j < 4; j++) {
1085 if (phi[j] < phimin) phimin = phi[j];
1086 if (phi[j] > phimax) phimax = phi[j];
1087 }
1088 }
1089 }
1090
1091 if (phimin > phimax) return fix_wcserr[status];
1092
1093 // Any changes needed?
1094 if (phimin >= -180.0 && phimax <= 180.0) return FIXERR_NO_CHANGE;
1095
1096
1097 // Compute the new reference pixel coordinates.
1098 double phi0 = (phimin + phimax) / 2.0;
1099 double theta0 = 0.0;
1100
1101 double x, y;
1102 if ((status = prjs2x(&(wcs->cel.prj), 1, 1, 1, 1, &phi0, &theta0, &x, &y,
1103 stat))) {
1104 if (status == PRJERR_BAD_PARAM) {
1105 status = FIXERR_BAD_PARAM;
1106 } else {
1107 status = FIXERR_NO_REF_PIX_COORD;
1108 }
1109 return wcserr_set(WCSFIX_ERRMSG(status));
1110 }
1111
1112 for (int k = 0; k < wcs->naxis; k++) {
1113 img[0][k] = 0.0;
1114 }
1115 img[0][wcs->lng] = x;
1116 img[0][wcs->lat] = y;
1117
1118 if ((status = linx2p(&(wcs->lin), 1, 0, img[0], pix[0]))) {
1119 return wcserr_set(WCSFIX_ERRMSG(fix_linerr[status]));
1120 }
1121
1122
1123 // Compute celestial coordinates at the new reference pixel.
1124 if ((status = wcsp2s(wcs, 1, 0, pix[0], img[0], phi, theta, world[0],
1125 stat))) {
1126 return fix_wcserr[status];
1127 }
1128
1129 // Compute native coordinates of the celestial pole.
1130 double lng = 0.0;
1131 double lat = 90.0;
1132 (void)sphs2x(wcs->cel.euler, 1, 1, 1, 1, &lng, &lat, phi, theta);
1133
1134 wcs->crpix[wcs->lng] = pix[0][wcs->lng];
1135 wcs->crpix[wcs->lat] = pix[0][wcs->lat];
1136 wcs->crval[wcs->lng] = world[0][wcs->lng];
1137 wcs->crval[wcs->lat] = world[0][wcs->lat];
1138 wcs->lonpole = phi[0] - phi0;
1139
1140 return wcsset(wcs);
1141 }
1142
1143 //----------------------------------------------------------------------------
1144
1145 // Helper function used only by wcspcx().
1146 static int unscramble(int n, int mapto[], int step, int type, void *vptr);
1147
wcspcx(struct wcsprm * wcs,int dopc,int permute,double rotn[2])1148 int wcspcx(
1149 struct wcsprm *wcs,
1150 int dopc,
1151 int permute,
1152 double rotn[2])
1153
1154 {
1155 static const char *function = "wcspcx";
1156
1157 // Initialize if required.
1158 if (wcs == 0x0) return FIXERR_NULL_POINTER;
1159 struct wcserr **err = &(wcs->err);
1160
1161 int status;
1162 if (wcs->flag != WCSSET) {
1163 if ((status = wcsset(wcs))) return fix_wcserr[status];
1164 }
1165
1166 // Check for CDi_j usage.
1167 double *wcscd = wcs->cd;
1168 if ((wcs->altlin & 1) || !(wcs->altlin & 2)) {
1169 if ((wcs->altlin & 1) && dopc == 1) {
1170 // Recompose PCi_j + CDELTi.
1171 wcscd = wcs->pc;
1172 } else {
1173 return wcserr_set(WCSERR_SET(FIXERR_BAD_PARAM),
1174 "CDi_j is not used in this coordinate representation");
1175 }
1176 }
1177
1178 // Check for sequent distortions.
1179 if (wcs->lin.disseq) {
1180 return wcserr_set(WCSERR_SET(FIXERR_BAD_COORD_TRANS),
1181 "Cannot handle coordinate descriptions containing sequent distortions");
1182 }
1183
1184
1185 // Allocate memory in bulk for two nxn matrices.
1186 int naxis = wcs->naxis;
1187 double *mem;
1188 if ((mem = calloc(2*naxis*naxis, sizeof(double))) == 0x0) {
1189 return wcserr_set(WCSFIX_ERRMSG(FIXERR_MEMORY));
1190 }
1191
1192 double *mat = mem;
1193 double *inv = mem + naxis*naxis;
1194
1195 // Construct the transpose of CDi_j with each element squared.
1196 double *matij = mat;
1197 for (int i = 0; i < naxis; i++) {
1198 double *cdji = wcscd + i;
1199 for (int j = 0; j < naxis; j++) {
1200 *(matij++) = (*cdji) * (*cdji);
1201 cdji += naxis;
1202 }
1203 }
1204
1205 // Invert the matrix.
1206 if ((status = matinv(naxis, mat, inv))) {
1207 return wcserr_set(WCSERR_SET(FIXERR_SINGULAR_MTX),
1208 "No solution for CDi_j matrix decomposition");
1209 }
1210
1211 // Apply scaling.
1212 double *invij = inv;
1213 double *pcij = wcs->pc;
1214 double *cdij = wcscd;
1215 for (int i = 0; i < naxis; i++) {
1216 double scl = 0.0;
1217 for (int j = 0; j < naxis; j++) {
1218 scl += *(invij++);
1219 }
1220
1221 scl = sqrt(scl);
1222 wcs->cdelt[i] /= scl;
1223
1224 for (int j = 0; j < naxis; j++) {
1225 *(pcij++) = *(cdij++) * scl;
1226 }
1227 }
1228
1229 // mapto[i] records where row i of PCi_j should move to.
1230 int *mapto = 0x0;
1231 if ((mapto = (int*)malloc(naxis*sizeof(int))) == 0x0) {
1232 free(mem);
1233 return wcserr_set(WCSFIX_ERRMSG(FIXERR_MEMORY));
1234 }
1235
1236 for (int i = 0; i < naxis; i++) {
1237 mapto[i] = -1;
1238 }
1239
1240 // Ensure that latitude always follows longitude.
1241 if (wcs->lng >= 0 && wcs->lat >= 0) {
1242 double *pci = wcs->pc + naxis*wcs->lng;
1243
1244 // Take the first non-zero element in the row.
1245 for (int j = 0; j < naxis; j++) {
1246 if (fabs(pci[j]) != 0.0) {
1247 mapto[wcs->lng] = j;
1248 break;
1249 }
1250 }
1251
1252 if (mapto[wcs->lng] == naxis-1) {
1253 mapto[wcs->lng]--;
1254 }
1255
1256 mapto[wcs->lat] = mapto[wcs->lng] + 1;
1257 }
1258
1259 // Fill in the rest of the row permutation map.
1260 for (int j = 0; j < naxis; j++) {
1261 // Column j.
1262 double *pcij = wcs->pc + j;
1263 double colmax = 0.0;
1264
1265 // Look down the column to find the absolute maximum element.
1266 for (int i = 0; i < naxis; i++, pcij += naxis) {
1267 if (!(mapto[i] < 0)) {
1268 // This row is already mapped.
1269 continue;
1270 }
1271
1272 if (fabs(*pcij) > colmax) {
1273 mapto[i] = j;
1274 colmax = fabs(*pcij);
1275 }
1276 }
1277 }
1278
1279
1280 // Fix the sign of CDELTi. Celestial axes are special, otherwise diagonal
1281 // elements of the correctly permuted matrix should be positive.
1282 for (int i = 0; i < naxis; i++) {
1283 int chsgn;
1284 double *pci = wcs->pc + naxis*i;
1285
1286 // Celestial axes are special.
1287 if (i == wcs->lng) {
1288 // Longitude axis - force CDELTi < 0.0.
1289 chsgn = (wcs->cdelt[i] > 0.0);
1290 } else if (i == wcs->lat) {
1291 // Latitude axis - force CDELTi > 0.0.
1292 chsgn = (wcs->cdelt[i] < 0.0);
1293 } else {
1294 chsgn = (pci[mapto[i]] < 0.0);
1295 }
1296
1297 if (chsgn) {
1298 wcs->cdelt[i] = -wcs->cdelt[i];
1299
1300 for (int j = 0; j < naxis; j++) {
1301 // Test needed to prevent negative zeros.
1302 if (pci[j] != 0.0) {
1303 pci[j] = -pci[j];
1304 }
1305 }
1306 }
1307 }
1308
1309 free(mem);
1310
1311 // Setting bit 3 in wcsprm::altlin stops wcsset() from reconstructing
1312 // PCi_j and CDELTi from CDi_j.
1313 wcs->altlin |= 8;
1314
1315
1316 // Compute rotation angle of each basis vector of the celestial axes.
1317 if (rotn) {
1318 if (wcs->lng < 0 || wcs->lat < 0) {
1319 // No celestial axes.
1320 rotn[0] = 0.0;
1321 rotn[1] = 0.0;
1322
1323 } else {
1324 double x, y;
1325 x = wcs->pc[naxis*wcs->lng + mapto[wcs->lng]];
1326 y = wcs->pc[naxis*wcs->lat + mapto[wcs->lng]];
1327 rotn[0] = atan2d(y, x);
1328
1329 y = -wcs->pc[naxis*wcs->lng + mapto[wcs->lat]];
1330 x = wcs->pc[naxis*wcs->lat + mapto[wcs->lat]];
1331 rotn[1] = atan2d(y, x);
1332 }
1333 }
1334
1335
1336 // Permute rows?
1337 if (permute) {
1338 // Check whether there's anything to unscramble.
1339 int scrambled = 0;
1340 for (int i = 0; i < naxis; i++) {
1341 if (mapto[i] != i) {
1342 scrambled = 1;
1343 break;
1344 }
1345 }
1346
1347 if (scrambled) {
1348 for (int i = 0; i < naxis; i++) {
1349 // Do columns of the PCi_ja matrix.
1350 if (unscramble(naxis, mapto, naxis, 1, wcs->pc + i)) goto cleanup;
1351 }
1352 if (unscramble(naxis, mapto, 1, 1, wcs->cdelt)) goto cleanup;
1353 if (unscramble(naxis, mapto, 1, 1, wcs->crval)) goto cleanup;
1354 if (unscramble(naxis, mapto, 1, 2, wcs->cunit)) goto cleanup;
1355 if (unscramble(naxis, mapto, 1, 2, wcs->ctype)) goto cleanup;
1356
1357 for (int ipv = 0; ipv < wcs->npv; ipv++) {
1358 // Noting that PVi_ma axis numbers are 1-relative.
1359 int i = wcs->pv[ipv].i - 1;
1360 wcs->pv[ipv].i = mapto[i] + 1;
1361 }
1362
1363 for (int ips = 0; ips < wcs->nps; ips++) {
1364 // Noting that PSi_ma axis numbers are 1-relative.
1365 int i = wcs->ps[ips].i - 1;
1366 wcs->ps[ips].i = mapto[i] + 1;
1367 }
1368
1369 if (wcs->altlin & 2) {
1370 for (int i = 0; i < naxis; i++) {
1371 // Do columns of the CDi_ja matrix.
1372 if (unscramble(naxis, mapto, naxis, 1, wcs->cd + i)) goto cleanup;
1373 }
1374 }
1375
1376 if (wcs->altlin & 4) {
1377 if (unscramble(naxis, mapto, 1, 1, wcs->crota)) goto cleanup;
1378 }
1379
1380 if (unscramble(naxis, mapto, 1, 3, wcs->colax)) goto cleanup;
1381 if (unscramble(naxis, mapto, 1, 2, wcs->cname)) goto cleanup;
1382 if (unscramble(naxis, mapto, 1, 1, wcs->crder)) goto cleanup;
1383 if (unscramble(naxis, mapto, 1, 1, wcs->csyer)) goto cleanup;
1384 if (unscramble(naxis, mapto, 1, 1, wcs->czphs)) goto cleanup;
1385 if (unscramble(naxis, mapto, 1, 1, wcs->cperi)) goto cleanup;
1386
1387 // Coordinate lookup tables.
1388 for (int itab = 0; itab < wcs->ntab; itab++) {
1389 for (int m = 0; m < wcs->tab[itab].M; m++) {
1390 int i = wcs->tab[itab].map[m];
1391 wcs->tab[itab].map[m] = mapto[i];
1392 }
1393 }
1394
1395 for (int iwtb = 0; iwtb < wcs->nwtb; iwtb++) {
1396 int i = wcs->wtb[iwtb].i;
1397 wcs->wtb[iwtb].i = mapto[i];
1398 }
1399
1400 // Distortions? No. Prior distortions operate on pixel coordinates and
1401 // therefore are not permuted, and sequent distortions are not handled.
1402 }
1403 }
1404
1405 free(mapto);
1406
1407 // Reset the struct.
1408 if ((status = wcsset(wcs))) return fix_wcserr[status];
1409
1410 return FIXERR_SUCCESS;
1411
1412 cleanup:
1413 if (mapto) free(mapto);
1414 return wcserr_set(WCSFIX_ERRMSG(FIXERR_MEMORY));
1415 }
1416
1417
unscramble(int n,int mapto[],int step,int type,void * vptr)1418 int unscramble(
1419 int n,
1420 int mapto[],
1421 int step,
1422 int type,
1423 void *vptr)
1424
1425 {
1426 if (step == 0) step = 1;
1427
1428 if (type == 1) {
1429 double *dval = (double *)vptr;
1430 double *dtmp;
1431 if ((dtmp = (double *)malloc(n*sizeof(double))) == 0x0) {
1432 return 1;
1433 }
1434
1435 for (int i = 0; i < n; i++) {
1436 dtmp[mapto[i]] = dval[i*step];
1437 }
1438
1439 for (int i = 0; i < n; i++) {
1440 dval[i*step] = dtmp[i];
1441 }
1442
1443 free(dtmp);
1444
1445 } else if (type == 2) {
1446 char (*cval)[72] = (char (*)[72])vptr;
1447 char (*ctmp)[72];
1448 if ((ctmp = (char (*)[72])malloc(n*72*sizeof(char))) == 0x0) {
1449 return 1;
1450 }
1451
1452 for (int i = 0; i < n; i++) {
1453 memcpy(ctmp[mapto[i]], cval[i], 72);
1454 }
1455
1456 for (int i = 0; i < n; i++) {
1457 memcpy(cval[i], ctmp[i], 72);
1458 }
1459
1460 free(ctmp);
1461
1462 } else if (type == 3) {
1463 int *ival = (int *)vptr;
1464 int *itmp;
1465 if ((itmp = (int *)malloc(n*sizeof(int))) == 0x0) {
1466 return 1;
1467 }
1468
1469 for (int i = 0; i < n; i++) {
1470 itmp[mapto[i]] = ival[i];
1471 }
1472
1473 for (int i = 0; i < n; i++) {
1474 ival[i] = itmp[i];
1475 }
1476
1477 free(itmp);
1478 }
1479
1480 return 0;
1481 }
1482