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: tab.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 "wcserr.h"
31 #include "wcsmath.h"
32 #include "wcsprintf.h"
33 #include "wcsutil.h"
34 #include "tab.h"
35
36 const int TABSET = 137;
37
38 // Map status return value to message.
39 const char *tab_errmsg[] = {
40 "Success",
41 "Null tabprm pointer passed",
42 "Memory allocation failed",
43 "Invalid tabular parameters",
44 "One or more of the x coordinates were invalid",
45 "One or more of the world coordinates were invalid"};
46
47 // Convenience macro for invoking wcserr_set().
48 #define TAB_ERRMSG(status) WCSERR_SET(status), tab_errmsg[status]
49
50 //----------------------------------------------------------------------------
51
tabini(int alloc,int M,const int K[],struct tabprm * tab)52 int tabini(int alloc, int M, const int K[], struct tabprm *tab)
53
54 {
55 static const char *function = "tabini";
56
57 if (tab == 0x0) return TABERR_NULL_POINTER;
58
59 // Initialize error message handling.
60 if (tab->flag == -1) {
61 tab->err = 0x0;
62 }
63 struct wcserr **err = &(tab->err);
64 wcserr_clear(err);
65
66
67 if (M <= 0) {
68 return wcserr_set(WCSERR_SET(TABERR_BAD_PARAMS),
69 "M must be positive, got %d", M);
70 }
71
72 // Determine the total number of elements in the coordinate array.
73 int N;
74 if (K) {
75 N = M;
76
77 for (int m = 0; m < M; m++) {
78 if (K[m] < 0) {
79 return wcserr_set(WCSERR_SET(TABERR_BAD_PARAMS),
80 "Invalid tabular parameters: Each element of K must be "
81 "non-negative, got %d", K[m]);
82 }
83
84 N *= K[m];
85 }
86
87 } else {
88 // Axis lengths as yet unknown.
89 N = 0;
90 }
91
92
93 // Initialize memory management.
94 if (tab->flag == -1 || tab->m_flag != TABSET) {
95 if (tab->flag == -1) {
96 tab->sense = 0x0;
97 tab->p0 = 0x0;
98 tab->delta = 0x0;
99 tab->extrema = 0x0;
100 tab->set_M = 0;
101 }
102
103 tab->m_flag = 0;
104 tab->m_M = 0;
105 tab->m_N = 0;
106 tab->m_K = 0x0;
107 tab->m_map = 0x0;
108 tab->m_crval = 0x0;
109 tab->m_index = 0x0;
110 tab->m_indxs = 0x0;
111 tab->m_coord = 0x0;
112
113 } else {
114 // Clear any outstanding signals set by wcstab().
115 for (int m = 0; m < tab->m_M; m++) {
116 if (tab->m_indxs[m] == (double *)0x1) tab->m_indxs[m] = 0x0;
117 }
118
119 if (tab->m_coord == (double *)0x1) tab->m_coord = 0x0;
120 }
121
122
123 // Allocate memory for arrays if required.
124 if (alloc ||
125 tab->K == 0x0 ||
126 tab->map == 0x0 ||
127 tab->crval == 0x0 ||
128 tab->index == 0x0 ||
129 tab->coord == 0x0) {
130
131 // Was sufficient allocated previously?
132 if (tab->m_flag == TABSET && (tab->m_M < M || tab->m_N < N)) {
133 // No, free it.
134 tabfree(tab);
135 }
136
137 if (alloc || tab->K == 0x0) {
138 if (tab->m_K) {
139 // In case the caller fiddled with it.
140 tab->K = tab->m_K;
141
142 } else {
143 if (!(tab->K = calloc(M, sizeof(int)))) {
144 return wcserr_set(TAB_ERRMSG(TABERR_MEMORY));
145 }
146
147 tab->m_flag = TABSET;
148 tab->m_M = M;
149 tab->m_K = tab->K;
150 }
151 }
152
153 if (alloc || tab->map == 0x0) {
154 if (tab->m_map) {
155 // In case the caller fiddled with it.
156 tab->map = tab->m_map;
157
158 } else {
159 if (!(tab->map = calloc(M, sizeof(int)))) {
160 return wcserr_set(TAB_ERRMSG(TABERR_MEMORY));
161 }
162
163 tab->m_flag = TABSET;
164 tab->m_M = M;
165 tab->m_map = tab->map;
166 }
167 }
168
169 if (alloc || tab->crval == 0x0) {
170 if (tab->m_crval) {
171 // In case the caller fiddled with it.
172 tab->crval = tab->m_crval;
173
174 } else {
175 if (!(tab->crval = calloc(M, sizeof(double)))) {
176 return wcserr_set(TAB_ERRMSG(TABERR_MEMORY));
177 }
178
179 tab->m_flag = TABSET;
180 tab->m_M = M;
181 tab->m_crval = tab->crval;
182 }
183 }
184
185 if (alloc || tab->index == 0x0) {
186 if (tab->m_index) {
187 // In case the caller fiddled with it.
188 tab->index = tab->m_index;
189
190 } else {
191 if (!(tab->index = calloc(M, sizeof(double *)))) {
192 return wcserr_set(TAB_ERRMSG(TABERR_MEMORY));
193 }
194
195 tab->m_flag = TABSET;
196 tab->m_M = M;
197 tab->m_N = N;
198 tab->m_index = tab->index;
199
200 if (!(tab->m_indxs = calloc(M, sizeof(double *)))) {
201 return wcserr_set(TAB_ERRMSG(TABERR_MEMORY));
202 }
203
204 // Recall that calloc() initializes these pointers to zero.
205 if (K) {
206 for (int m = 0; m < M; m++) {
207 if (K[m]) {
208 if (!(tab->index[m] = calloc(K[m], sizeof(double)))) {
209 return wcserr_set(TAB_ERRMSG(TABERR_MEMORY));
210 }
211
212 tab->m_indxs[m] = tab->index[m];
213 }
214 }
215 }
216 }
217 }
218
219 if (alloc || tab->coord == 0x0) {
220 if (tab->m_coord) {
221 // In case the caller fiddled with it.
222 tab->coord = tab->m_coord;
223
224 } else if (N) {
225 if (!(tab->coord = calloc(N, sizeof(double)))) {
226 return wcserr_set(TAB_ERRMSG(TABERR_MEMORY));
227 }
228
229 tab->m_flag = TABSET;
230 tab->m_M = M;
231 tab->m_N = N;
232 tab->m_coord = tab->coord;
233 }
234 }
235 }
236
237 tab->flag = 0;
238 tab->M = M;
239
240 // Set defaults.
241 for (int m = 0; m < M; m++) {
242 tab->map[m] = -1;
243 tab->crval[m] = 0.0;
244
245 if (K) {
246 tab->K[m] = K[m];
247 double *dp;
248 if ((dp = tab->index[m])) {
249 for (int k = 0; k < K[m]; k++) {
250 *(dp++) = k;
251 }
252 }
253 } else {
254 tab->K[m] = 0;
255 }
256 }
257
258 // Initialize the coordinate array.
259 for (double *dp = tab->coord; dp < tab->coord + N; dp++) {
260 *dp = UNDEFINED;
261 }
262
263 return 0;
264 }
265
266 //----------------------------------------------------------------------------
267
tabmem(struct tabprm * tab)268 int tabmem(struct tabprm *tab)
269
270 {
271 static const char *function = "tabmem";
272
273 if (tab == 0x0) return TABERR_NULL_POINTER;
274 struct wcserr **err = &(tab->err);
275
276 if (tab->M == 0 || tab->K == 0x0) {
277 // Should have been set by this time.
278 return wcserr_set(WCSERR_SET(TABERR_MEMORY),
279 "Null pointers in tabprm struct");
280 }
281
282
283 int M = tab->M;
284 int N = tab->M;
285 for (int m = 0; m < M; m++) {
286 if (tab->K[m] < 0) {
287 return wcserr_set(WCSERR_SET(TABERR_BAD_PARAMS),
288 "Invalid tabular parameters: Each element of K must be "
289 "non-negative, got %d", M);
290 }
291
292 N *= tab->K[m];
293 }
294
295
296 if (tab->m_M == 0) {
297 tab->m_M = M;
298 } else if (tab->m_M < M) {
299 // Only possible if the user changed M.
300 return wcserr_set(WCSERR_SET(TABERR_MEMORY),
301 "tabprm struct inconsistent");
302 }
303
304 if (tab->m_N == 0) {
305 tab->m_N = N;
306 } else if (tab->m_N < N) {
307 // Only possible if the user changed K[].
308 return wcserr_set(WCSERR_SET(TABERR_MEMORY),
309 "tabprm struct inconsistent");
310 }
311
312 if (tab->m_K == 0x0) {
313 if ((tab->m_K = tab->K)) {
314 tab->m_flag = TABSET;
315 }
316 }
317
318 if (tab->m_map == 0x0) {
319 if ((tab->m_map = tab->map)) {
320 tab->m_flag = TABSET;
321 }
322 }
323
324 if (tab->m_crval == 0x0) {
325 if ((tab->m_crval = tab->crval)) {
326 tab->m_flag = TABSET;
327 }
328 }
329
330 if (tab->m_index == 0x0) {
331 if ((tab->m_index = tab->index)) {
332 tab->m_flag = TABSET;
333 }
334 }
335
336 for (int m = 0; m < tab->m_M; m++) {
337 if (tab->m_indxs[m] == 0x0 || tab->m_indxs[m] == (double *)0x1) {
338 if ((tab->m_indxs[m] = tab->index[m])) {
339 tab->m_flag = TABSET;
340 }
341 }
342 }
343
344 if (tab->m_coord == 0x0 || tab->m_coord == (double *)0x1) {
345 if ((tab->m_coord = tab->coord)) {
346 tab->m_flag = TABSET;
347 }
348 }
349
350 tab->flag = 0;
351
352 return 0;
353 }
354
355 //----------------------------------------------------------------------------
356
tabcpy(int alloc,const struct tabprm * tabsrc,struct tabprm * tabdst)357 int tabcpy(int alloc, const struct tabprm *tabsrc, struct tabprm *tabdst)
358
359 {
360 static const char *function = "tabcpy";
361
362 int status;
363
364 if (tabsrc == 0x0) return TABERR_NULL_POINTER;
365 if (tabdst == 0x0) return TABERR_NULL_POINTER;
366 struct wcserr **err = &(tabdst->err);
367
368 int M = tabsrc->M;
369 if (M <= 0) {
370 return wcserr_set(WCSERR_SET(TABERR_BAD_PARAMS),
371 "M must be positive, got %d", M);
372 }
373
374 if ((status = tabini(alloc, M, tabsrc->K, tabdst))) {
375 return status;
376 }
377
378 int N = M;
379 for (int m = 0; m < M; m++) {
380 tabdst->map[m] = tabsrc->map[m];
381 tabdst->crval[m] = tabsrc->crval[m];
382 N *= tabsrc->K[m];
383 }
384
385 double *dstp, *srcp;
386 for (int m = 0; m < M; m++) {
387 if ((srcp = tabsrc->index[m])) {
388 dstp = tabdst->index[m];
389 for (int k = 0; k < tabsrc->K[m]; k++) {
390 *(dstp++) = *(srcp++);
391 }
392 }
393 }
394
395 srcp = tabsrc->coord;
396 dstp = tabdst->coord;
397 for (int n = 0; n < N; n++) {
398 *(dstp++) = *(srcp++);
399 }
400
401 return 0;
402 }
403
404 //----------------------------------------------------------------------------
405
tabcmp(int dummy,double tol,const struct tabprm * tab1,const struct tabprm * tab2,int * equal)406 int tabcmp(
407 int dummy,
408 double tol,
409 const struct tabprm *tab1,
410 const struct tabprm *tab2,
411 int *equal)
412
413 {
414 // Avert nuisance compiler warnings about unused parameters.
415 (void)dummy;
416
417 if (tab1 == 0x0) return TABERR_NULL_POINTER;
418 if (tab2 == 0x0) return TABERR_NULL_POINTER;
419 if (equal == 0x0) return TABERR_NULL_POINTER;
420
421 *equal = 0;
422
423 if (tab1->M != tab2->M) {
424 return 0;
425 }
426
427 int M = tab1->M;
428
429 if (!wcsutil_intEq(M, tab1->K, tab2->K) ||
430 !wcsutil_intEq(M, tab1->map, tab2->map) ||
431 !wcsutil_dblEq(M, tol, tab1->crval, tab2->crval)) {
432 return 0;
433 }
434
435 int N = M;
436 for (int m = 0; m < M; m++) {
437 if (!wcsutil_dblEq(tab1->K[m], tol, tab1->index[m], tab2->index[m])) {
438 return 0;
439 }
440
441 N *= tab1->K[m];
442 }
443
444 if (!wcsutil_dblEq(N, tol, tab1->coord, tab2->coord)) {
445 return 0;
446 }
447
448 *equal = 1;
449
450 return 0;
451 }
452
453
454 //----------------------------------------------------------------------------
455
tabfree(struct tabprm * tab)456 int tabfree(struct tabprm *tab)
457
458 {
459 if (tab == 0x0) return TABERR_NULL_POINTER;
460
461 if (tab->flag != -1) {
462 // Clear any outstanding signals set by wcstab().
463 for (int m = 0; m < tab->m_M; m++) {
464 if (tab->m_indxs[m] == (double *)0x1) tab->m_indxs[m] = 0x0;
465 }
466
467 if (tab->m_coord == (double *)0x1) tab->m_coord = 0x0;
468
469 // Free memory allocated by tabini().
470 if (tab->m_flag == TABSET) {
471 if (tab->K == tab->m_K) tab->K = 0x0;
472 if (tab->map == tab->m_map) tab->map = 0x0;
473 if (tab->crval == tab->m_crval) tab->crval = 0x0;
474 if (tab->index == tab->m_index) tab->index = 0x0;
475 if (tab->coord == tab->m_coord) tab->coord = 0x0;
476
477 if (tab->m_K) free(tab->m_K);
478 if (tab->m_map) free(tab->m_map);
479 if (tab->m_crval) free(tab->m_crval);
480
481 if (tab->m_index) {
482 for (int m = 0; m < tab->m_M; m++) {
483 if (tab->m_indxs[m]) free(tab->m_indxs[m]);
484 }
485 free(tab->m_index);
486 free(tab->m_indxs);
487 }
488
489 if (tab->m_coord) free(tab->m_coord);
490 }
491
492 // Free memory allocated by tabset().
493 if (tab->sense) free(tab->sense);
494 if (tab->p0) free(tab->p0);
495 if (tab->delta) free(tab->delta);
496 if (tab->extrema) free(tab->extrema);
497 }
498
499 tab->m_flag = 0;
500 tab->m_M = 0;
501 tab->m_N = 0;
502 tab->m_K = 0x0;
503 tab->m_map = 0x0;
504 tab->m_crval = 0x0;
505 tab->m_index = 0x0;
506 tab->m_indxs = 0x0;
507 tab->m_coord = 0x0;
508
509 tab->sense = 0x0;
510 tab->p0 = 0x0;
511 tab->delta = 0x0;
512 tab->extrema = 0x0;
513 tab->set_M = 0;
514
515 wcserr_clear(&(tab->err));
516
517 tab->flag = 0;
518
519 return 0;
520 }
521
522 //----------------------------------------------------------------------------
523
tabsize(const struct tabprm * tab,int sizes[2])524 int tabsize(const struct tabprm *tab, int sizes[2])
525
526 {
527 if (tab == 0x0) {
528 sizes[0] = sizes[1] = 0;
529 return TABERR_SUCCESS;
530 }
531
532 // Base size, in bytes.
533 sizes[0] = sizeof(struct tabprm);
534
535 // Total size of allocated memory, in bytes.
536 sizes[1] = 0;
537
538 int exsizes[2];
539 int M = tab->M;
540
541 // tabprm::K[];
542 sizes[1] += M * sizeof(int);
543
544 // tabprm::map[];
545 sizes[1] += M * sizeof(int);
546
547 // tabprm::crval[];
548 sizes[1] += M * sizeof(double);
549
550 // tabprm::index[] and tabprm::m_indxs;
551 sizes[1] += 2*M * sizeof(double *);
552 for (int m = 0; m < M; m++) {
553 if (tab->index[m]) {
554 sizes[1] += tab->K[m] * sizeof(double);
555 }
556 }
557
558 // tabprm::coord[];
559 sizes[1] += M * tab->nc * sizeof(double);
560
561 // tab::err[].
562 wcserr_size(tab->err, exsizes);
563 sizes[1] += exsizes[0] + exsizes[1];
564
565 // The remaining arrays are allocated by tabset().
566 if (tab->flag != TABSET) {
567 return TABERR_SUCCESS;
568 }
569
570 // tabprm::sense[].
571 if (tab->sense) {
572 sizes[1] += M * sizeof(int);
573 }
574
575 // tabprm::p0[].
576 if (tab->p0) {
577 sizes[1] += M * sizeof(int);
578 }
579
580 // tabprm::delta[].
581 if (tab->delta) {
582 sizes[1] += M * sizeof(double);
583 }
584
585 // tabprm::extrema[].
586 int ne = (tab->nc / tab->K[0]) * 2 * M;
587 sizes[1] += ne * sizeof(double);
588
589 return TABERR_SUCCESS;
590 }
591
592 //----------------------------------------------------------------------------
593
tabprt(const struct tabprm * tab)594 int tabprt(const struct tabprm *tab)
595
596 {
597 char *cp, text[128];
598 double *dp;
599
600 if (tab == 0x0) return TABERR_NULL_POINTER;
601
602 if (tab->flag != TABSET) {
603 wcsprintf("The tabprm struct is UNINITIALIZED.\n");
604 return 0;
605 }
606
607 wcsprintf(" flag: %d\n", tab->flag);
608 wcsprintf(" M: %d\n", tab->M);
609
610 // Array dimensions.
611 WCSPRINTF_PTR(" K: ", tab->K, "\n");
612 wcsprintf(" ");
613 for (int m = 0; m < tab->M; m++) {
614 wcsprintf("%6d", tab->K[m]);
615 }
616 wcsprintf("\n");
617
618 // Map vector.
619 WCSPRINTF_PTR(" map: ", tab->map, "\n");
620 wcsprintf(" ");
621 for (int m = 0; m < tab->M; m++) {
622 wcsprintf("%6d", tab->map[m]);
623 }
624 wcsprintf("\n");
625
626 // Reference index value.
627 WCSPRINTF_PTR(" crval: ", tab->crval, "\n");
628 wcsprintf(" ");
629 for (int m = 0; m < tab->M; m++) {
630 wcsprintf(" %#- 11.5g", tab->crval[m]);
631 }
632 wcsprintf("\n");
633
634 // Index vectors.
635 WCSPRINTF_PTR(" index: ", tab->index, "\n");
636 for (int m = 0; m < tab->M; m++) {
637 wcsprintf(" index[%d]: ", m);
638 WCSPRINTF_PTR("", tab->index[m], "");
639 if (tab->index[m]) {
640 for (int k = 0; k < tab->K[m]; k++) {
641 if (k%5 == 0) {
642 wcsprintf("\n ");
643 }
644 wcsprintf(" %#- 11.5g", tab->index[m][k]);
645 }
646 wcsprintf("\n");
647 }
648 }
649
650 // Coordinate array.
651 WCSPRINTF_PTR(" coord: ", tab->coord, "\n");
652 dp = tab->coord;
653 for (int n = 0; n < tab->nc; n++) {
654 // Array index.
655 int j = n;
656 cp = text;
657 for (int m = 0; m < tab->M; m++) {
658 int nd = (tab->K[m] < 10) ? 1 : 2;
659 sprintf(cp, ",%*d", nd, j % tab->K[m] + 1);
660 j /= tab->K[m];
661 cp += strlen(cp);
662 }
663
664 wcsprintf(" (*%s)", text);
665 for (int m = 0; m < tab->M; m++) {
666 wcsprintf(" %#- 11.5g", *(dp++));
667 }
668 wcsprintf("\n");
669 }
670
671 wcsprintf(" nc: %d\n", tab->nc);
672
673 WCSPRINTF_PTR(" sense: ", tab->sense, "\n");
674 if (tab->sense) {
675 wcsprintf(" ");
676 for (int m = 0; m < tab->M; m++) {
677 wcsprintf("%6d", tab->sense[m]);
678 }
679 wcsprintf("\n");
680 }
681
682 WCSPRINTF_PTR(" p0: ", tab->p0, "\n");
683 if (tab->p0) {
684 wcsprintf(" ");
685 for (int m = 0; m < tab->M; m++) {
686 wcsprintf("%6d", tab->p0[m]);
687 }
688 wcsprintf("\n");
689 }
690
691 WCSPRINTF_PTR(" delta: ", tab->delta, "\n");
692 if (tab->delta) {
693 wcsprintf(" ");
694 for (int m = 0; m < tab->M; m++) {
695 wcsprintf(" %#- 11.5g", tab->delta[m]);
696 }
697 wcsprintf("\n");
698 }
699
700 WCSPRINTF_PTR(" extrema: ", tab->extrema, "\n");
701 dp = tab->extrema;
702 for (int n = 0; n < tab->nc/tab->K[0]; n++) {
703 // Array index.
704 int j = n;
705 cp = text;
706 *cp = '\0';
707 for (int m = 1; m < tab->M; m++) {
708 int nd = (tab->K[m] < 10) ? 1 : 2;
709 sprintf(cp, ",%*d", nd, j % tab->K[m] + 1);
710 j /= tab->K[m];
711 cp += strlen(cp);
712 }
713
714 wcsprintf(" (*,*%s)", text);
715 for (int m = 0; m < 2*tab->M; m++) {
716 if (m == tab->M) wcsprintf("-> ");
717 wcsprintf(" %#- 11.5g", *(dp++));
718 }
719 wcsprintf("\n");
720 }
721
722 WCSPRINTF_PTR(" err: ", tab->err, "\n");
723 if (tab->err) {
724 wcserr_prt(tab->err, " ");
725 }
726
727 // Memory management.
728 wcsprintf(" m_flag: %d\n", tab->m_flag);
729 wcsprintf(" m_M: %d\n", tab->m_M);
730 wcsprintf(" m_N: %d\n", tab->m_N);
731
732 WCSPRINTF_PTR(" m_K: ", tab->m_K, "");
733 if (tab->m_K == tab->K) wcsprintf(" (= K)");
734 wcsprintf("\n");
735
736 WCSPRINTF_PTR(" m_map: ", tab->m_map, "");
737 if (tab->m_map == tab->map) wcsprintf(" (= map)");
738 wcsprintf("\n");
739
740 WCSPRINTF_PTR(" m_crval: ", tab->m_crval, "");
741 if (tab->m_crval == tab->crval) wcsprintf(" (= crval)");
742 wcsprintf("\n");
743
744 WCSPRINTF_PTR(" m_index: ", tab->m_index, "");
745 if (tab->m_index == tab->index) wcsprintf(" (= index)");
746 wcsprintf("\n");
747 for (int m = 0; m < tab->M; m++) {
748 wcsprintf(" m_indxs[%d]: ", m);
749 WCSPRINTF_PTR("", tab->m_indxs[m], "");
750 if (tab->m_indxs[m] == tab->index[m]) wcsprintf(" (= index[%d])", m);
751 wcsprintf("\n");
752 }
753
754 WCSPRINTF_PTR(" m_coord: ", tab->m_coord, "");
755 if (tab->m_coord == tab->coord) wcsprintf(" (= coord)");
756 wcsprintf("\n");
757
758 return 0;
759 }
760
761 //----------------------------------------------------------------------------
762
tabperr(const struct tabprm * tab,const char * prefix)763 int tabperr(const struct tabprm *tab, const char *prefix)
764
765 {
766 if (tab == 0x0) return TABERR_NULL_POINTER;
767
768 if (tab->err) {
769 wcserr_prt(tab->err, prefix);
770 }
771
772 return 0;
773 }
774
775 //----------------------------------------------------------------------------
776
tabset(struct tabprm * tab)777 int tabset(struct tabprm *tab)
778
779 {
780 static const char *function = "tabset";
781
782 if (tab == 0x0) return TABERR_NULL_POINTER;
783 struct wcserr **err = &(tab->err);
784
785 // Check the number of tabular coordinate axes.
786 int M = tab->M;
787 if (M < 1) {
788 return wcserr_set(WCSERR_SET(TABERR_BAD_PARAMS),
789 "Invalid tabular parameters: M must be positive, got %d", M);
790 }
791
792 // Check the axis lengths.
793 if (!tab->K) {
794 return wcserr_set(WCSERR_SET(TABERR_MEMORY),
795 "Null pointers in tabprm struct");
796 }
797
798 tab->nc = 1;
799 for (int m = 0; m < M; m++) {
800 if (tab->K[m] < 1) {
801 return wcserr_set(WCSERR_SET(TABERR_BAD_PARAMS),
802 "Invalid tabular parameters: Each element of K must be positive, "
803 "got %d", tab->K[m]);
804 }
805
806 // Number of coordinate vectors in the coordinate array.
807 tab->nc *= tab->K[m];
808 }
809
810 // Check that the map vector is sensible.
811 if (!tab->map) {
812 return wcserr_set(WCSERR_SET(TABERR_MEMORY),
813 "Null pointers in tabprm struct");
814 }
815
816 for (int m = 0; m < M; m++) {
817 int i = tab->map[m];
818 if (i < 0) {
819 return wcserr_set(WCSERR_SET(TABERR_BAD_PARAMS),
820 "Invalid tabular parameters: Each element of map must be "
821 "non-negative, got %d", i);
822 }
823 }
824
825 // Check memory allocation for the remaining vectors.
826 if (!tab->crval || !tab->index || !tab->coord) {
827 return wcserr_set(WCSERR_SET(TABERR_MEMORY),
828 "Null pointers in tabprm struct");
829 }
830
831 // Take memory if signalled to by wcstab().
832 for (int m = 0; m < tab->m_M; m++) {
833 if (tab->m_indxs[m] == (double *)0x1 &&
834 (tab->m_indxs[m] = tab->index[m])) {
835 tab->m_flag = TABSET;
836 }
837 }
838
839 if (tab->m_coord == (double *)0x1 &&
840 (tab->m_coord = tab->coord)) {
841 tab->m_flag = TABSET;
842 }
843
844
845 // Allocate memory for work vectors.
846 if (tab->flag != TABSET || tab->set_M < M) {
847 // Free memory that may have been allocated previously.
848 if (tab->sense) free(tab->sense);
849 if (tab->p0) free(tab->p0);
850 if (tab->delta) free(tab->delta);
851 if (tab->extrema) free(tab->extrema);
852
853 // Allocate memory for internal arrays.
854 if (!(tab->sense = calloc(M, sizeof(int)))) {
855 return wcserr_set(TAB_ERRMSG(TABERR_MEMORY));
856 }
857
858 if (!(tab->p0 = calloc(M, sizeof(int)))) {
859 free(tab->sense);
860 return wcserr_set(TAB_ERRMSG(TABERR_MEMORY));
861 }
862
863 if (!(tab->delta = calloc(M, sizeof(double)))) {
864 free(tab->sense);
865 free(tab->p0);
866 return wcserr_set(TAB_ERRMSG(TABERR_MEMORY));
867 }
868
869 int ne = (tab->nc / tab->K[0]) * 2 * M;
870 if (!(tab->extrema = calloc(ne, sizeof(double)))) {
871 free(tab->sense);
872 free(tab->p0);
873 free(tab->delta);
874 return wcserr_set(TAB_ERRMSG(TABERR_MEMORY));
875 }
876
877 tab->set_M = M;
878 }
879
880 // Check that the index vectors are monotonic.
881 int *Km = tab->K;
882 for (int m = 0; m < M; m++, Km++) {
883 tab->sense[m] = 0;
884
885 if (*Km > 1) {
886 double *Psi;
887 if ((Psi = tab->index[m]) == 0x0) {
888 // Default indexing.
889 tab->sense[m] = 1;
890
891 } else {
892 for (int k = 0; k < *Km-1; k++) {
893 switch (tab->sense[m]) {
894 case 0:
895 if (Psi[k] < Psi[k+1]) {
896 // Monotonic increasing.
897 tab->sense[m] = 1;
898 } else if (Psi[k] > Psi[k+1]) {
899 // Monotonic decreasing.
900 tab->sense[m] = -1;
901 }
902 break;
903
904 case 1:
905 if (Psi[k] > Psi[k+1]) {
906 // Should be monotonic increasing.
907 free(tab->sense);
908 free(tab->p0);
909 free(tab->delta);
910 free(tab->extrema);
911 return wcserr_set(WCSERR_SET(TABERR_BAD_PARAMS),
912 "Invalid tabular parameters: Index vectors are not "
913 "monotonically increasing");
914 }
915 break;
916
917 case -1:
918 if (Psi[k] < Psi[k+1]) {
919 // Should be monotonic decreasing.
920 free(tab->sense);
921 free(tab->p0);
922 free(tab->delta);
923 free(tab->extrema);
924 return wcserr_set(WCSERR_SET(TABERR_BAD_PARAMS),
925 "Invalid tabular parameters: Index vectors are not "
926 "monotonically decreasing");
927 }
928 break;
929 }
930 }
931 }
932
933 if (tab->sense[m] == 0) {
934 free(tab->sense);
935 free(tab->p0);
936 free(tab->delta);
937 free(tab->extrema);
938 return wcserr_set(WCSERR_SET(TABERR_BAD_PARAMS),
939 "Invalid tabular parameters: Index vectors are not monotonic");
940 }
941 }
942 }
943
944 // Find the extremal values of the coordinate elements in each row.
945 double *dcrd = tab->coord;
946 double *dmin = tab->extrema;
947 double *dmax = tab->extrema + M;
948 for (int ic = 0; ic < tab->nc; ic += tab->K[0]) {
949 for (int m = 0; m < M; m++, dcrd++) {
950 if (tab->K[0] > 1) {
951 // Extrapolate a little before the start of the row.
952 double dPsi;
953 double *Psi = tab->index[0];
954 if (Psi == 0x0) {
955 dPsi = 1.0;
956 } else {
957 dPsi = Psi[1] - Psi[0];
958 }
959
960 double dval = *dcrd;
961 if (dPsi != 0.0) {
962 dval -= 0.5 * (*(dcrd+M) - *dcrd)/dPsi;
963 }
964
965 *(dmax+m) = *(dmin+m) = dval;
966 } else {
967 *(dmax+m) = *(dmin+m) = *dcrd;
968 }
969 }
970
971 dcrd -= M;
972 for (int i = 0; i < tab->K[0]; i++) {
973 for (int m = 0; m < M; m++, dcrd++) {
974 if (*(dmax+m) < *dcrd) *(dmax+m) = *dcrd;
975 if (*(dmin+m) > *dcrd) *(dmin+m) = *dcrd;
976
977 if (tab->K[0] > 1 && i == tab->K[0]-1) {
978 // Extrapolate a little beyond the end of the row.
979 double dPsi;
980 double *Psi = tab->index[0];
981 if (Psi == 0x0) {
982 dPsi = 1.0;
983 } else {
984 dPsi = Psi[i] - Psi[i-1];
985 }
986
987 double dval = *dcrd;
988 if (dPsi != 0.0) {
989 dval += 0.5 * (*dcrd - *(dcrd-M))/dPsi;
990 }
991
992 if (*(dmax+m) < dval) *(dmax+m) = dval;
993 if (*(dmin+m) > dval) *(dmin+m) = dval;
994 }
995 }
996 }
997
998 dmin += 2*M;
999 dmax += 2*M;
1000 }
1001
1002 tab->flag = TABSET;
1003
1004 return 0;
1005 }
1006
1007 //----------------------------------------------------------------------------
1008
tabx2s(struct tabprm * tab,int ncoord,int nelem,const double x[],double world[],int stat[])1009 int tabx2s(
1010 struct tabprm *tab,
1011 int ncoord,
1012 int nelem,
1013 const double x[],
1014 double world[],
1015 int stat[])
1016
1017 {
1018 static const char *function = "tabx2s";
1019
1020 int status;
1021
1022 if (tab == 0x0) return TABERR_NULL_POINTER;
1023 struct wcserr **err = &(tab->err);
1024
1025 // Initialize if required.
1026 if (tab->flag != TABSET) {
1027 if ((status = tabset(tab))) return status;
1028 }
1029
1030 // This is used a lot.
1031 int M = tab->M;
1032
1033 status = 0;
1034 register const double *xp = x;
1035 register double *wp = world;
1036 register int *statp = stat;
1037 for (int n = 0; n < ncoord; n++) {
1038 // Determine the indexes.
1039 int *Km = tab->K;
1040 for (int m = 0; m < M; m++, Km++) {
1041 // N.B. psi_m and Upsilon_m are 1-relative FITS indexes.
1042 int i = tab->map[m];
1043 double psi_m = *(xp+i) + tab->crval[m];
1044
1045 double *Psi = tab->index[m];
1046 double upsilon;
1047 if (Psi == 0x0) {
1048 // Default indexing is simple.
1049 upsilon = psi_m;
1050
1051 } else {
1052 // To ease confusion, decrement Psi so that we can use 1-relative
1053 // C array indexing to match the 1-relative FITS indexing.
1054 Psi--;
1055
1056 if (*Km == 1) {
1057 // Index vector is degenerate.
1058 if (Psi[1]-0.5 <= psi_m && psi_m <= Psi[1]+0.5) {
1059 upsilon = psi_m;
1060 } else {
1061 *statp = 1;
1062 status = wcserr_set(TAB_ERRMSG(TABERR_BAD_X));
1063 goto next;
1064 }
1065
1066 } else {
1067 // Interpolate in the indexing vector.
1068 int k;
1069 if (tab->sense[m] == 1) {
1070 // Monotonic increasing index values.
1071 if (psi_m < Psi[1]) {
1072 if (Psi[1] - 0.5*(Psi[2]-Psi[1]) <= psi_m) {
1073 // Allow minor extrapolation.
1074 k = 1;
1075
1076 } else {
1077 // Index is out of range.
1078 *statp = 1;
1079 status = wcserr_set(TAB_ERRMSG(TABERR_BAD_X));
1080 goto next;
1081 }
1082
1083 } else if (Psi[*Km] < psi_m) {
1084 if (psi_m <= Psi[*Km] + 0.5*(Psi[*Km]-Psi[*Km-1])) {
1085 // Allow minor extrapolation.
1086 k = *Km - 1;
1087
1088 } else {
1089 // Index is out of range.
1090 *statp = 1;
1091 status = wcserr_set(TAB_ERRMSG(TABERR_BAD_X));
1092 goto next;
1093 }
1094
1095 } else {
1096 for (k = 1; k < *Km; k++) {
1097 if (psi_m < Psi[k]) {
1098 continue;
1099 }
1100 if (Psi[k] == psi_m && psi_m < Psi[k+1]) {
1101 break;
1102 }
1103 if (Psi[k] < psi_m && psi_m <= Psi[k+1]) {
1104 break;
1105 }
1106 }
1107 }
1108
1109 } else {
1110 // Monotonic decreasing index values.
1111 if (psi_m > Psi[1]) {
1112 if (Psi[1] + 0.5*(Psi[1]-Psi[2]) >= psi_m) {
1113 // Allow minor extrapolation.
1114 k = 1;
1115
1116 } else {
1117 // Index is out of range.
1118 *statp = 1;
1119 status = wcserr_set(TAB_ERRMSG(TABERR_BAD_X));
1120 goto next;
1121 }
1122
1123 } else if (psi_m < Psi[*Km]) {
1124 if (Psi[*Km] - 0.5*(Psi[*Km-1]-Psi[*Km]) <= psi_m) {
1125 // Allow minor extrapolation.
1126 k = *Km - 1;
1127
1128 } else {
1129 // Index is out of range.
1130 *statp = 1;
1131 status = wcserr_set(TAB_ERRMSG(TABERR_BAD_X));
1132 goto next;
1133 }
1134
1135 } else {
1136 for (k = 1; k < *Km; k++) {
1137 if (psi_m > Psi[k]) {
1138 continue;
1139 }
1140 if (Psi[k] == psi_m && psi_m > Psi[k+1]) {
1141 break;
1142 }
1143 if (Psi[k] > psi_m && psi_m >= Psi[k+1]) {
1144 break;
1145 }
1146 }
1147 }
1148 }
1149
1150 upsilon = k + (psi_m - Psi[k]) / (Psi[k+1] - Psi[k]);
1151 }
1152 }
1153
1154 if (upsilon < 0.5 || upsilon > *Km + 0.5) {
1155 // Index out of range.
1156 *statp = 1;
1157 status = wcserr_set(TAB_ERRMSG(TABERR_BAD_X));
1158 goto next;
1159 }
1160
1161 // Fiducial array indices and fractional offset.
1162 // p1 is 1-relative while tab::p0 is 0-relative.
1163 int p1 = (int)floor(upsilon);
1164 tab->p0[m] = p1 - 1;
1165 tab->delta[m] = upsilon - p1;
1166
1167 if (p1 == 0) {
1168 // Extrapolation below p1 == 1.
1169 tab->p0[m] += 1;
1170 tab->delta[m] -= 1.0;
1171 } else if (p1 == *Km && *Km > 1) {
1172 // Extrapolation above p1 == K_m.
1173 tab->p0[m] -= 1;
1174 tab->delta[m] += 1.0;
1175 }
1176 }
1177
1178
1179 // Now interpolate in the coordinate array; the M-dimensional linear
1180 // interpolation algorithm is described in Sect. 3.4 of WCS Paper IV.
1181 for (int m = 0; m < M; m++) {
1182 int i = tab->map[m];
1183 *(wp+i) = 0.0;
1184 }
1185
1186 // Loop over the 2^M vertices surrounding P.
1187 int nv = 1 << M;
1188 for (int iv = 0; iv < nv; iv++) {
1189 // Locate vertex in the coordinate array and compute its weight.
1190 int offset = 0;
1191 double wgt = 1.0;
1192 for (int m = M-1; m >= 0; m--) {
1193 offset *= tab->K[m];
1194 offset += tab->p0[m];
1195 if (iv & (1 << m)) {
1196 if (tab->K[m] > 1) offset++;
1197 wgt *= tab->delta[m];
1198 } else {
1199 wgt *= 1.0 - tab->delta[m];
1200 }
1201 }
1202
1203 if (wgt == 0.0) continue;
1204
1205 // Add the contribution from this vertex to each element.
1206 double *coord = tab->coord + offset*M;
1207 for (int m = 0; m < M; m++) {
1208 int i = tab->map[m];
1209 *(wp+i) += *(coord++) * wgt;
1210 }
1211
1212 if (wgt == 1.0) break;
1213 }
1214
1215 *statp = 0;
1216
1217 next:
1218 xp += nelem;
1219 wp += nelem;
1220 statp++;
1221 }
1222
1223 return status;
1224 }
1225
1226 //----------------------------------------------------------------------------
1227
1228 // Helper functions used only by tabs2x().
1229 static int tabedge(struct tabprm *);
1230 static int tabrow(struct tabprm *, const double *);
1231 static int tabvox(struct tabprm *, const double *, int, double **,
1232 unsigned int *);
1233
tabs2x(struct tabprm * tab,int ncoord,int nelem,const double world[],double x[],int stat[])1234 int tabs2x(
1235 struct tabprm* tab,
1236 int ncoord,
1237 int nelem,
1238 const double world[],
1239 double x[],
1240 int stat[])
1241
1242 {
1243 static const char *function = "tabs2x";
1244
1245 int status;
1246
1247 if (tab == 0x0) return TABERR_NULL_POINTER;
1248 struct wcserr **err = &(tab->err);
1249
1250 // Initialize if required.
1251 if (tab->flag != TABSET) {
1252 if ((status = tabset(tab))) return status;
1253 }
1254
1255 // This is used a lot.
1256 int M = tab->M;
1257
1258 double **tabcoord = 0x0;
1259 int nv = 0;
1260 if (M > 1) {
1261 nv = 1 << M;
1262 tabcoord = calloc(nv, sizeof(double *));
1263 }
1264
1265
1266 status = 0;
1267 register const double *wp = world;
1268 register double *xp = x;
1269 register int *statp = stat;
1270 for (int n = 0; n < ncoord; n++) {
1271 // Locate this coordinate in the coordinate array.
1272 int edge = 0;
1273 for (int m = 0; m < M; m++) {
1274 tab->p0[m] = 0;
1275 }
1276
1277 int ic;
1278 for (ic = 0; ic < tab->nc; ic++) {
1279 if (tab->p0[0] == 0) {
1280 // New row, could it contain a solution?
1281 if (edge || tabrow(tab, wp)) {
1282 // No, skip it.
1283 ic += tab->K[0];
1284 if (1 < M) {
1285 tab->p0[1]++;
1286 edge = tabedge(tab);
1287 }
1288
1289 // Because ic will be incremented when the loop is reentered.
1290 ic--;
1291 continue;
1292 }
1293 }
1294
1295 if (M == 1) {
1296 // Deal with the one-dimensional case separately for efficiency.
1297 double w = wp[tab->map[0]];
1298 if (w == tab->coord[0]) {
1299 tab->p0[0] = 0;
1300 tab->delta[0] = 0.0;
1301 break;
1302
1303 } else if (ic < tab->nc - 1) {
1304 if (((tab->coord[ic] <= w && w <= tab->coord[ic+1]) ||
1305 (tab->coord[ic] >= w && w >= tab->coord[ic+1])) &&
1306 (tab->index[0] == 0x0 ||
1307 tab->index[0][ic] != tab->index[0][ic+1])) {
1308 tab->p0[0] = ic;
1309 tab->delta[0] = (w - tab->coord[ic]) /
1310 (tab->coord[ic+1] - tab->coord[ic]);
1311 break;
1312 }
1313 }
1314
1315 } else {
1316 // Multi-dimensional tables are harder.
1317 if (!edge) {
1318 // Addresses of the coordinates for each corner of the "voxel".
1319 for (int iv = 0; iv < nv; iv++) {
1320 int offset = 0;
1321 for (int m = M-1; m >= 0; m--) {
1322 offset *= tab->K[m];
1323 offset += tab->p0[m];
1324 if ((iv & (1 << m)) && (tab->K[m] > 1)) offset++;
1325 }
1326 tabcoord[iv] = tab->coord + offset*M;
1327 }
1328
1329 if (tabvox(tab, wp, 0, tabcoord, 0x0) == 0) {
1330 // Found a solution.
1331 break;
1332 }
1333 }
1334
1335 // Next voxel.
1336 tab->p0[0]++;
1337 edge = tabedge(tab);
1338 }
1339 }
1340
1341
1342 if (ic == tab->nc) {
1343 // Coordinate not found; allow minor extrapolation.
1344 if (M == 1) {
1345 // Should there be a solution?
1346 double w = wp[tab->map[0]];
1347 if (tab->extrema[0] <= w && w <= tab->extrema[1]) {
1348 double *dcrd = tab->coord;
1349 for (int i = 0; i < 2; i++) {
1350 if (i) dcrd += tab->K[0] - 2;
1351
1352 double delta = (w - *dcrd) / (*(dcrd+1) - *dcrd);
1353
1354 if (i == 0) {
1355 if (-0.5 <= delta && delta <= 0.0) {
1356 tab->p0[0] = 0;
1357 tab->delta[0] = delta;
1358 ic = 0;
1359 break;
1360 }
1361 } else {
1362 if (1.0 <= delta && delta <= 1.5) {
1363 tab->p0[0] = tab->K[0] - 1;
1364 tab->delta[0] = delta - 1.0;
1365 ic = 0;
1366 }
1367 }
1368 }
1369 }
1370
1371 } else {
1372 // Multi-dimensional tables.
1373 // >>> TBD <<<
1374 }
1375 }
1376
1377
1378 if (ic == tab->nc) {
1379 // Coordinate not found.
1380 *statp = 1;
1381 status = wcserr_set(TAB_ERRMSG(TABERR_BAD_WORLD));
1382
1383 } else {
1384 // Determine the intermediate world coordinates.
1385 int *Km = tab->K;
1386 for (int m = 0; m < M; m++, Km++) {
1387 // N.B. Upsilon_m and psi_m are 1-relative FITS indexes.
1388 double upsilon = (tab->p0[m] + 1) + tab->delta[m];
1389
1390 if (upsilon < 0.5 || upsilon > *Km + 0.5) {
1391 // Index out of range.
1392 *statp = 1;
1393 status = wcserr_set(TAB_ERRMSG(TABERR_BAD_WORLD));
1394
1395 } else {
1396 // Do inverse lookup of the index vector.
1397 double *Psi = tab->index[m];
1398 double psi_m;
1399 if (Psi == 0x0) {
1400 // Default indexing.
1401 psi_m = upsilon;
1402
1403 } else {
1404 // Decrement Psi and use 1-relative C array indexing to match the
1405 // 1-relative FITS indexing.
1406 Psi--;
1407
1408 if (*Km == 1) {
1409 // Degenerate index vector.
1410 psi_m = Psi[1];
1411 } else {
1412 int k = (int)(upsilon);
1413 psi_m = Psi[k];
1414 if (k < *Km) {
1415 psi_m += (upsilon - k) * (Psi[k+1] - Psi[k]);
1416 }
1417 }
1418 }
1419
1420 xp[tab->map[m]] = psi_m - tab->crval[m];
1421 }
1422 }
1423 *statp = 0;
1424 }
1425
1426 wp += nelem;
1427 xp += nelem;
1428 statp++;
1429 }
1430
1431 if (tabcoord) free(tabcoord);
1432
1433 return status;
1434 }
1435
1436 /*----------------------------------------------------------------------------
1437 * Convenience routine to check whether tabprm::p0 has been incremented beyond
1438 * the end of an index vector and if so move it to the start of the next one.
1439 * Returns 1 if tabprm::p0 is sitting at the end of any non-degenerate index
1440 * vector.
1441 *---------------------------------------------------------------------------*/
1442
tabedge(struct tabprm * tab)1443 int tabedge(struct tabprm* tab)
1444
1445 {
1446 int edge = 0;
1447
1448 for (int m = 0; m < tab->M; m++) {
1449 if (tab->p0[m] == tab->K[m]) {
1450 // p0 has been incremented beyond the end of an index vector, point it
1451 // to the next one.
1452 tab->p0[m] = 0;
1453 if (m < tab->M-1) {
1454 tab->p0[m+1]++;
1455 }
1456 } else if (tab->p0[m] == tab->K[m]-1 && tab->K[m] > 1) {
1457 // p0 is sitting at the end of a non-degenerate index vector.
1458 edge = 1;
1459 }
1460 }
1461
1462 return edge;
1463 }
1464
1465 /*----------------------------------------------------------------------------
1466 * Quick test to see whether the world coordinate indicated by wp could lie
1467 * somewhere along (or near) the row of the image indexed by tabprm::p0.
1468 * Return 0 if so, 1 otherwise.
1469 *
1470 * tabprm::p0 selects a particular row of the image, p0[0] being ignored (i.e.
1471 * treated as zero). Adjacent rows that delimit a row of "voxels" are formed
1472 * by incrementing elements other than p0[0] in all binary combinations. N.B.
1473 * these are not the same as the voxels (pixels) that are indexed by, and
1474 * centred on, integral pixel coordinates in FITS.
1475 *
1476 * To see why it is necessary to examine the adjacent rows, consider the 2-D
1477 * case where the first world coordinate element is constant along each row.
1478 * If the first element of wp has value 0.5, and its value in the row indexed
1479 * by p0 has value 0, and in the next row it has value 1, then it is clear that
1480 * the solution lies in neither row but somewhere between them. Thus both rows
1481 * will be involved in finding the solution.
1482 *
1483 * tabprm::extrema is the address of the first element of a 1-D array that
1484 * records the minimum and maximum value of each element of the coordinate
1485 * vector in each row of the coordinate array, treated as though it were
1486 * defined as
1487 *
1488 * double extrema[K_M]...[K_2][2][M]
1489 *
1490 * The minimum is recorded in the first element of the compressed K_1
1491 * dimension, then the maximum.
1492 *---------------------------------------------------------------------------*/
1493
tabrow(struct tabprm * tab,const double * wp)1494 int tabrow(struct tabprm* tab, const double *wp)
1495
1496 {
1497 const double tol = 1e-10;
1498
1499 int M = tab->M;
1500
1501 // The number of corners in a "voxel". We need examine only half this
1502 // number of rows. The extra factor of two will be used to select between
1503 // the minimal and maximal values in each row.
1504 unsigned int nv = 1 << M;
1505
1506 unsigned int eq = 0;
1507 unsigned int lt = 0;
1508 unsigned int gt = 0;
1509 for (unsigned int iv = 0; iv < nv; iv++) {
1510 // Find the index into tabprm::extrema for this row.
1511 int offset = 0;
1512 for (int m = M-1; m > 0; m--) {
1513 offset *= tab->K[m];
1514 offset += tab->p0[m];
1515
1516 // Select the row.
1517 if (iv & (1 << m)) {
1518 if (tab->K[m] > 1) offset++;
1519 }
1520 }
1521
1522 // The K_1 dimension has length 2 (see prologue).
1523 offset *= 2;
1524
1525 // Select the minimum on even numbered iterations, else the maximum.
1526 if (iv & 1) offset++;
1527
1528 // The last dimension has length M (see prologue).
1529 offset *= M;
1530
1531 // Address of the extremal elements (min or max) for this row.
1532 double *cp = tab->extrema + offset;
1533
1534 // For each coordinate element, we only need to find one row where its
1535 // minimum value is less than that of wp, and one row where the maximum
1536 // value is greater. That doesn't mean that there is a solution, only
1537 // that there might be.
1538 for (int m = 0; m < M; m++, cp++) {
1539 // Apply the axis mapping.
1540 double w = wp[tab->map[m]];
1541
1542 // Finally the test itself; set bits in the bitmask.
1543 if (fabs(*cp - w) < tol) {
1544 eq |= (1 << m);
1545 } else if (*cp < w) {
1546 lt |= (1 << m);
1547 } else if (*cp > w) {
1548 gt |= (1 << m);
1549 }
1550 }
1551
1552 // Have all bits been switched on?
1553 if ((lt | eq) == nv-1 && (gt | eq) == nv-1) {
1554 // A solution could lie within this row of voxels.
1555 return 0;
1556 }
1557 }
1558
1559 // No solution in this row.
1560 return 1;
1561 }
1562
1563 /*----------------------------------------------------------------------------
1564 * Does the world coordinate indicated by wp lie within the voxel indexed by
1565 * tabprm::p0? If so, do a binary chop of the interior of the voxel to find
1566 * it and return 0, with tabprm::delta set to the solution. Else return 1.
1567 *
1568 * As in tabrow(), a "voxel" is formed by incrementing the elements of
1569 * tabprm::p0 in all binary combinations. Note that these are not the same as
1570 * the voxels (pixels) that are indexed by, and centred on, integral pixel
1571 * coordinates in FITS.
1572 *
1573 * tabvox() calls itself recursively. When called from outside, level, being
1574 * the level of recursion, should be given as zero. tabcoord is an array
1575 * holding the addresses of the coordinates for each corner of the voxel.
1576 * vox is the address of a work array (vox2) used during recursive calls to
1577 * dissect the voxel. It is ignored when tabvox() is called from outside
1578 * (level == 0).
1579 *
1580 * It is assumed that the image dimensions are no greater than 32.
1581 ----------------------------------------------------------------------------*/
1582
tabvox(struct tabprm * tab,const double * wp,int level,double ** tabcoord,unsigned int * vox)1583 int tabvox(
1584 struct tabprm* tab,
1585 const double *wp,
1586 int level,
1587 double **tabcoord,
1588 unsigned int *vox)
1589
1590 {
1591 const double tol = 1e-10;
1592
1593 int M = tab->M;
1594
1595 // The number of corners in a voxel.
1596 unsigned int nv = 1 << M;
1597
1598 double dv = 1.0;
1599 for (int i = 0; i < level; i++) {
1600 dv /= 2.0;
1601 }
1602
1603 // Could the coordinate lie within this voxel (level == 0) or sub-voxel
1604 // (level > 0)? We use the fact that with linear interpolation the
1605 // coordinate elements are extremal in a corner and test each one.
1606 unsigned int lt = 0;
1607 unsigned int gt = 0;
1608 unsigned int eq = 0;
1609 for (unsigned int iv = 0; iv < nv; iv++) {
1610 // Select a corner of the sub-voxel.
1611 double coord[32];
1612 for (int m = 0; m < M; m++) {
1613 coord[m] = 0.0;
1614 tab->delta[m] = level ? dv*vox[m] : 0.0;
1615
1616 if (iv & (1 << m)) {
1617 tab->delta[m] += dv;
1618 }
1619 }
1620
1621 // Compute the coordinates of this corner of the sub-voxel by linear
1622 // interpolation using the weighting algorithm described in Sect. 3.4 of
1623 // WCS Paper IV.
1624 for (unsigned int jv = 0; jv < nv; jv++) {
1625 // Find the weight for this corner of the parent voxel.
1626 double wgt = 1.0;
1627 for (int m = 0; m < M; m++) {
1628 if (jv & (1 << m)) {
1629 wgt *= tab->delta[m];
1630 } else {
1631 wgt *= 1.0 - tab->delta[m];
1632 }
1633 }
1634
1635 if (wgt == 0.0) continue;
1636
1637 // Add its contribution to each coordinate element.
1638 double *cp = tabcoord[jv];
1639 for (int m = 0; m < M; m++) {
1640 coord[m] += *(cp++) * wgt;
1641 }
1642
1643 if (wgt == 1.0) break;
1644 }
1645
1646 // Coordinate elements are minimal or maximal in a corner.
1647 unsigned int et = 0;
1648 for (int m = 0; m < M; m++) {
1649 // Apply the axis mapping.
1650 double w = wp[tab->map[m]];
1651
1652 // Finally the test itself; set bits in the bitmask.
1653 if (fabs(coord[m] - w) < tol) {
1654 et |= (1 << m);
1655 } else if (coord[m] < w) {
1656 lt |= (1 << m);
1657 } else if (coord[m] > w) {
1658 gt |= (1 << m);
1659 }
1660 }
1661
1662 if (et == nv-1) {
1663 // We've stumbled across a solution in this corner of the sub-voxel.
1664 return 0;
1665 }
1666
1667 eq |= et;
1668 }
1669
1670 // Could the coordinate lie within this sub-voxel?
1671 if ((lt | eq) == nv-1 && (gt | eq) == nv-1) {
1672 // Yes it could, but does it?
1673
1674 // Is it time to stop the recursion?
1675 if (level == 31) {
1676 // We have a solution, squeeze out the last bit of juice.
1677 dv /= 2.0;
1678 for (int m = 0; m < M; m++) {
1679 tab->delta[m] = dv * (2.0*vox[m] + 1.0);
1680 }
1681
1682 return 0;
1683 }
1684
1685 // Subdivide the sub-voxel and try again for each subdivision.
1686 for (unsigned int iv = 0; iv < nv; iv++) {
1687 // Select the subdivision.
1688 unsigned int vox2[32];
1689 for (int m = 0; m < M; m++) {
1690 vox2[m] = level ? 2*vox[m] : 0;
1691 if (iv & (1 << m)) {
1692 vox2[m]++;
1693 }
1694 }
1695
1696 // Recurse.
1697 if (tabvox(tab, wp, level+1, tabcoord, vox2) == 0) {
1698 return 0;
1699 }
1700 }
1701 }
1702
1703 // No solution in this sub-voxel.
1704 return 1;
1705 }
1706