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