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