1 /*============================================================================ 2 WCSLIB 7.7 - an implementation of the FITS WCS standard. 3 Copyright (C) 1995-2021, Mark Calabretta 4 5 This file is part of WCSLIB. 6 7 WCSLIB is free software: you can redistribute it and/or modify it under the 8 terms of the GNU Lesser General Public License as published by the Free 9 Software Foundation, either version 3 of the License, or (at your option) 10 any later version. 11 12 WCSLIB is distributed in the hope that it will be useful, but WITHOUT ANY 13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for 15 more details. 16 17 You should have received a copy of the GNU Lesser General Public License 18 along with WCSLIB. If not, see http://www.gnu.org/licenses. 19 20 Author: Mark Calabretta, Australia Telescope National Facility, CSIRO. 21 http://www.atnf.csiro.au/people/Mark.Calabretta 22 $Id: tab.h,v 7.7 2021/07/12 06:36:49 mcalabre Exp $ 23 *============================================================================= 24 * 25 * WCSLIB 7.7 - C routines that implement the FITS World Coordinate System 26 * (WCS) standard. Refer to the README file provided with WCSLIB for an 27 * overview of the library. 28 * 29 * 30 * Summary of the tab routines 31 * --------------------------- 32 * Routines in this suite implement the part of the FITS World Coordinate 33 * System (WCS) standard that deals with tabular coordinates, i.e. coordinates 34 * that are defined via a lookup table, as described in 35 * 36 = "Representations of world coordinates in FITS", 37 = Greisen, E.W., & Calabretta, M.R. 2002, A&A, 395, 1061 (WCS Paper I) 38 = 39 = "Representations of spectral coordinates in FITS", 40 = Greisen, E.W., Calabretta, M.R., Valdes, F.G., & Allen, S.L. 41 = 2006, A&A, 446, 747 (WCS Paper III) 42 * 43 * These routines define methods to be used for computing tabular world 44 * coordinates from intermediate world coordinates (a linear transformation 45 * of image pixel coordinates), and vice versa. They are based on the tabprm 46 * struct which contains all information needed for the computations. The 47 * struct contains some members that must be set by the user, and others that 48 * are maintained by these routines, somewhat like a C++ class but with no 49 * encapsulation. 50 * 51 * tabini(), tabmem(), tabcpy(), and tabfree() are provided to manage the 52 * tabprm struct, tabsize() computes its total size including allocated memory, 53 * and tabprt() prints its contents. 54 * 55 * tabperr() prints the error message(s) (if any) stored in a tabprm struct. 56 * 57 * A setup routine, tabset(), computes intermediate values in the tabprm struct 58 * from parameters in it that were supplied by the user. The struct always 59 * needs to be set up by tabset() but it need not be called explicitly - refer 60 * to the explanation of tabprm::flag. 61 * 62 * tabx2s() and tabs2x() implement the WCS tabular coordinate transformations. 63 * 64 * Accuracy: 65 * --------- 66 * No warranty is given for the accuracy of these routines (refer to the 67 * copyright notice); intending users must satisfy for themselves their 68 * adequacy for the intended purpose. However, closure effectively to within 69 * double precision rounding error was demonstrated by test routine ttab.c 70 * which accompanies this software. 71 * 72 * 73 * tabini() - Default constructor for the tabprm struct 74 * ---------------------------------------------------- 75 * tabini() allocates memory for arrays in a tabprm struct and sets all members 76 * of the struct to default values. 77 * 78 * PLEASE NOTE: every tabprm struct should be initialized by tabini(), possibly 79 * repeatedly. On the first invokation, and only the first invokation, the 80 * flag member of the tabprm struct must be set to -1 to initialize memory 81 * management, regardless of whether tabini() will actually be used to allocate 82 * memory. 83 * 84 * Given: 85 * alloc int If true, allocate memory unconditionally for arrays in 86 * the tabprm struct. 87 * 88 * If false, it is assumed that pointers to these arrays 89 * have been set by the user except if they are null 90 * pointers in which case memory will be allocated for 91 * them regardless. (In other words, setting alloc true 92 * saves having to initalize these pointers to zero.) 93 * 94 * M int The number of tabular coordinate axes. 95 * 96 * K const int[] 97 * Vector of length M whose elements (K_1, K_2,... K_M) 98 * record the lengths of the axes of the coordinate array 99 * and of each indexing vector. M and K[] are used to 100 * determine the length of the various tabprm arrays and 101 * therefore the amount of memory to allocate for them. 102 * Their values are copied into the tabprm struct. 103 * 104 * It is permissible to set K (i.e. the address of the 105 * array) to zero which has the same effect as setting 106 * each element of K[] to zero. In this case no memory 107 * will be allocated for the index vectors or coordinate 108 * array in the tabprm struct. These together with the 109 * K vector must be set separately before calling 110 * tabset(). 111 * 112 * Given and returned: 113 * tab struct tabprm* 114 * Tabular transformation parameters. Note that, in 115 * order to initialize memory management tabprm::flag 116 * should be set to -1 when tab is initialized for the 117 * first time (memory leaks may result if it had already 118 * been initialized). 119 * 120 * Function return value: 121 * int Status return value: 122 * 0: Success. 123 * 1: Null tabprm pointer passed. 124 * 2: Memory allocation failed. 125 * 3: Invalid tabular parameters. 126 * 127 * For returns > 1, a detailed error message is set in 128 * tabprm::err if enabled, see wcserr_enable(). 129 * 130 * 131 * tabmem() - Acquire tabular memory 132 * --------------------------------- 133 * tabmem() takes control of memory allocated by the user for arrays in the 134 * tabprm struct. 135 * 136 * Given and returned: 137 * tab struct tabprm* 138 * Tabular transformation parameters. 139 * 140 * Function return value: 141 * int Status return value: 142 * 0: Success. 143 * 1: Null tabprm pointer passed. 144 * 2: Memory allocation failed. 145 * 146 * For returns > 1, a detailed error message is set in 147 * tabprm::err if enabled, see wcserr_enable(). 148 * 149 * 150 * tabcpy() - Copy routine for the tabprm struct 151 * --------------------------------------------- 152 * tabcpy() does a deep copy of one tabprm struct to another, using tabini() to 153 * allocate memory for its arrays if required. Only the "information to be 154 * provided" part of the struct is copied; a call to tabset() is required to 155 * set up the remainder. 156 * 157 * Given: 158 * alloc int If true, allocate memory unconditionally for arrays in 159 * the tabprm struct. 160 * 161 * If false, it is assumed that pointers to these arrays 162 * have been set by the user except if they are null 163 * pointers in which case memory will be allocated for 164 * them regardless. (In other words, setting alloc true 165 * saves having to initalize these pointers to zero.) 166 * 167 * tabsrc const struct tabprm* 168 * Struct to copy from. 169 * 170 * Given and returned: 171 * tabdst struct tabprm* 172 * Struct to copy to. tabprm::flag should be set to -1 173 * if tabdst was not previously initialized (memory leaks 174 * may result if it was previously initialized). 175 * 176 * Function return value: 177 * int Status return value: 178 * 0: Success. 179 * 1: Null tabprm pointer passed. 180 * 2: Memory allocation failed. 181 * 182 * For returns > 1, a detailed error message is set in 183 * tabprm::err (associated with tabdst) if enabled, see 184 * wcserr_enable(). 185 * 186 * 187 * tabcmp() - Compare two tabprm structs for equality 188 * -------------------------------------------------- 189 * tabcmp() compares two tabprm structs for equality. 190 * 191 * Given: 192 * cmp int A bit field controlling the strictness of the 193 * comparison. At present, this value must always be 0, 194 * indicating a strict comparison. In the future, other 195 * options may be added. 196 * 197 * tol double Tolerance for comparison of floating-point values. 198 * For example, for tol == 1e-6, all floating-point 199 * values in the structs must be equal to the first 6 200 * decimal places. A value of 0 implies exact equality. 201 * 202 * tab1 const struct tabprm* 203 * The first tabprm struct to compare. 204 * 205 * tab2 const struct tabprm* 206 * The second tabprm struct to compare. 207 * 208 * Returned: 209 * equal int* Non-zero when the given structs are equal. 210 * 211 * Function return value: 212 * int Status return value: 213 * 0: Success. 214 * 1: Null pointer passed. 215 * 216 * 217 * tabfree() - Destructor for the tabprm struct 218 * -------------------------------------------- 219 * tabfree() frees memory allocated for the tabprm arrays by tabini(). 220 * tabini() records the memory it allocates and tabfree() will only attempt to 221 * free this. 222 * 223 * PLEASE NOTE: tabfree() must not be invoked on a tabprm struct that was not 224 * initialized by tabini(). 225 * 226 * Returned: 227 * tab struct tabprm* 228 * Coordinate transformation parameters. 229 * 230 * Function return value: 231 * int Status return value: 232 * 0: Success. 233 * 1: Null tabprm pointer passed. 234 * 235 * 236 * tabsize() - Compute the size of a tabprm struct 237 * ----------------------------------------------- 238 * tabsize() computes the full size of a tabprm struct, including allocated 239 * memory. 240 * 241 * Given: 242 * tab const struct tabprm* 243 * Tabular transformation parameters. 244 * 245 * If NULL, the base size of the struct and the allocated 246 * size are both set to zero. 247 * 248 * Returned: 249 * sizes int[2] The first element is the base size of the struct as 250 * returned by sizeof(struct tabprm). The second element 251 * is the total allocated size, in bytes, assuming that 252 * the allocation was done by tabini(). This figure 253 * includes memory allocated for the constituent struct, 254 * tabprm::err. 255 * 256 * It is not an error for the struct not to have been set 257 * up via tabset(), which normally results in additional 258 * memory allocation. 259 * 260 * Function return value: 261 * int Status return value: 262 * 0: Success. 263 * 264 * 265 * tabprt() - Print routine for the tabprm struct 266 * ---------------------------------------------- 267 * tabprt() prints the contents of a tabprm struct using wcsprintf(). Mainly 268 * intended for diagnostic purposes. 269 * 270 * Given: 271 * tab const struct tabprm* 272 * Tabular transformation parameters. 273 * 274 * Function return value: 275 * int Status return value: 276 * 0: Success. 277 * 1: Null tabprm pointer passed. 278 * 279 * 280 * tabperr() - Print error messages from a tabprm struct 281 * ----------------------------------------------------- 282 * tabperr() prints the error message(s) (if any) stored in a tabprm struct. 283 * If there are no errors then nothing is printed. It uses wcserr_prt(), q.v. 284 * 285 * Given: 286 * tab const struct tabprm* 287 * Tabular transformation parameters. 288 * 289 * prefix const char * 290 * If non-NULL, each output line will be prefixed with 291 * this string. 292 * 293 * Function return value: 294 * int Status return value: 295 * 0: Success. 296 * 1: Null tabprm pointer passed. 297 * 298 * 299 * tabset() - Setup routine for the tabprm struct 300 * ----------------------------------------------- 301 * tabset() allocates memory for work arrays in the tabprm struct and sets up 302 * the struct according to information supplied within it. 303 * 304 * Note that this routine need not be called directly; it will be invoked by 305 * tabx2s() and tabs2x() if tabprm::flag is anything other than a predefined 306 * magic value. 307 * 308 * Given and returned: 309 * tab struct tabprm* 310 * Tabular transformation parameters. 311 * 312 * Function return value: 313 * int Status return value: 314 * 0: Success. 315 * 1: Null tabprm pointer passed. 316 * 3: Invalid tabular parameters. 317 * 318 * For returns > 1, a detailed error message is set in 319 * tabprm::err if enabled, see wcserr_enable(). 320 * 321 * 322 * tabx2s() - Pixel-to-world transformation 323 * ---------------------------------------- 324 * tabx2s() transforms intermediate world coordinates to world coordinates 325 * using coordinate lookup. 326 * 327 * Given and returned: 328 * tab struct tabprm* 329 * Tabular transformation parameters. 330 * 331 * Given: 332 * ncoord, 333 * nelem int The number of coordinates, each of vector length 334 * nelem. 335 * 336 * x const double[ncoord][nelem] 337 * Array of intermediate world coordinates, SI units. 338 * 339 * Returned: 340 * world double[ncoord][nelem] 341 * Array of world coordinates, in SI units. 342 * 343 * stat int[ncoord] 344 * Status return value status for each coordinate: 345 * 0: Success. 346 * 1: Invalid intermediate world coordinate. 347 * 348 * Function return value: 349 * int Status return value: 350 * 0: Success. 351 * 1: Null tabprm pointer passed. 352 * 3: Invalid tabular parameters. 353 * 4: One or more of the x coordinates were invalid, 354 * as indicated by the stat vector. 355 * 356 * For returns > 1, a detailed error message is set in 357 * tabprm::err if enabled, see wcserr_enable(). 358 * 359 * 360 * tabs2x() - World-to-pixel transformation 361 * ---------------------------------------- 362 * tabs2x() transforms world coordinates to intermediate world coordinates. 363 * 364 * Given and returned: 365 * tab struct tabprm* 366 * Tabular transformation parameters. 367 * 368 * Given: 369 * ncoord, 370 * nelem int The number of coordinates, each of vector length 371 * nelem. 372 * world const double[ncoord][nelem] 373 * Array of world coordinates, in SI units. 374 * 375 * Returned: 376 * x double[ncoord][nelem] 377 * Array of intermediate world coordinates, SI units. 378 * stat int[ncoord] 379 * Status return value status for each vector element: 380 * 0: Success. 381 * 1: Invalid world coordinate. 382 * 383 * Function return value: 384 * int Status return value: 385 * 0: Success. 386 * 1: Null tabprm pointer passed. 387 * 3: Invalid tabular parameters. 388 * 5: One or more of the world coordinates were 389 * invalid, as indicated by the stat vector. 390 * 391 * For returns > 1, a detailed error message is set in 392 * tabprm::err if enabled, see wcserr_enable(). 393 * 394 * 395 * tabprm struct - Tabular transformation parameters 396 * ------------------------------------------------- 397 * The tabprm struct contains information required to transform tabular 398 * coordinates. It consists of certain members that must be set by the user 399 * ("given") and others that are set by the WCSLIB routines ("returned"). Some 400 * of the latter are supplied for informational purposes while others are for 401 * internal use only. 402 * 403 * int flag 404 * (Given and returned) This flag must be set to zero whenever any of the 405 * following tabprm structure members are set or changed: 406 * 407 * - tabprm::M (q.v., not normally set by the user), 408 * - tabprm::K (q.v., not normally set by the user), 409 * - tabprm::map, 410 * - tabprm::crval, 411 * - tabprm::index, 412 * - tabprm::coord. 413 * 414 * This signals the initialization routine, tabset(), to recompute the 415 * returned members of the tabprm struct. tabset() will reset flag to 416 * indicate that this has been done. 417 * 418 * PLEASE NOTE: flag should be set to -1 when tabini() is called for the 419 * first time for a particular tabprm struct in order to initialize memory 420 * management. It must ONLY be used on the first initialization otherwise 421 * memory leaks may result. 422 * 423 * int M 424 * (Given or returned) Number of tabular coordinate axes. 425 * 426 * If tabini() is used to initialize the tabprm struct (as would normally 427 * be the case) then it will set M from the value passed to it as a 428 * function argument. The user should not subsequently modify it. 429 * 430 * int *K 431 * (Given or returned) Pointer to the first element of a vector of length 432 * tabprm::M whose elements (K_1, K_2,... K_M) record the lengths of the 433 * axes of the coordinate array and of each indexing vector. 434 * 435 * If tabini() is used to initialize the tabprm struct (as would normally 436 * be the case) then it will set K from the array passed to it as a 437 * function argument. The user should not subsequently modify it. 438 * 439 * int *map 440 * (Given) Pointer to the first element of a vector of length tabprm::M 441 * that defines the association between axis m in the M-dimensional 442 * coordinate array (1 <= m <= M) and the indices of the intermediate world 443 * coordinate and world coordinate arrays, x[] and world[], in the argument 444 * lists for tabx2s() and tabs2x(). 445 * 446 * When x[] and world[] contain the full complement of coordinate elements 447 * in image-order, as will usually be the case, then map[m-1] == i-1 for 448 * axis i in the N-dimensional image (1 <= i <= N). In terms of the FITS 449 * keywords 450 * 451 * map[PVi_3a - 1] == i - 1. 452 * 453 * However, a different association may result if x[], for example, only 454 * contains a (relevant) subset of intermediate world coordinate elements. 455 * For example, if M == 1 for an image with N > 1, it is possible to fill 456 * x[] with the relevant coordinate element with nelem set to 1. In this 457 * case map[0] = 0 regardless of the value of i. 458 * 459 * double *crval 460 * (Given) Pointer to the first element of a vector of length tabprm::M 461 * whose elements contain the index value for the reference pixel for each 462 * of the tabular coordinate axes. 463 * 464 * double **index 465 * (Given) Pointer to the first element of a vector of length tabprm::M of 466 * pointers to vectors of lengths (K_1, K_2,... K_M) of 0-relative indexes 467 * (see tabprm::K). 468 * 469 * The address of any or all of these index vectors may be set to zero, 470 * i.e. 471 * 472 = index[m] == 0; 473 * 474 * this is interpreted as default indexing, i.e. 475 * 476 = index[m][k] = k; 477 * 478 * double *coord 479 * (Given) Pointer to the first element of the tabular coordinate array, 480 * treated as though it were defined as 481 * 482 = double coord[K_M]...[K_2][K_1][M]; 483 * 484 * (see tabprm::K) i.e. with the M dimension varying fastest so that the 485 * M elements of a coordinate vector are stored contiguously in memory. 486 * 487 * int nc 488 * (Returned) Total number of coordinate vectors in the coordinate array 489 * being the product K_1 * K_2 * ... * K_M (see tabprm::K). 490 * 491 * int padding 492 * (An unused variable inserted for alignment purposes only.) 493 * 494 * int *sense 495 * (Returned) Pointer to the first element of a vector of length tabprm::M 496 * whose elements indicate whether the corresponding indexing vector is 497 * monotonic increasing (+1), or decreasing (-1). 498 * 499 * int *p0 500 * (Returned) Pointer to the first element of a vector of length tabprm::M 501 * of interpolated indices into the coordinate array such that Upsilon_m, 502 * as defined in Paper III, is equal to (p0[m] + 1) + tabprm::delta[m]. 503 * 504 * double *delta 505 * (Returned) Pointer to the first element of a vector of length tabprm::M 506 * of interpolated indices into the coordinate array such that Upsilon_m, 507 * as defined in Paper III, is equal to (tabprm::p0[m] + 1) + delta[m]. 508 * 509 * double *extrema 510 * (Returned) Pointer to the first element of an array that records the 511 * minimum and maximum value of each element of the coordinate vector in 512 * each row of the coordinate array, treated as though it were defined as 513 * 514 = double extrema[K_M]...[K_2][2][M] 515 * 516 * (see tabprm::K). The minimum is recorded in the first element of the 517 * compressed K_1 dimension, then the maximum. This array is used by the 518 * inverse table lookup function, tabs2x(), to speed up table searches. 519 * 520 * struct wcserr *err 521 * (Returned) If enabled, when an error status is returned, this struct 522 * contains detailed information about the error, see wcserr_enable(). 523 * 524 * int m_flag 525 * (For internal use only.) 526 * int m_M 527 * (For internal use only.) 528 * int m_N 529 * (For internal use only.) 530 * int set_M 531 * (For internal use only.) 532 * int m_K 533 * (For internal use only.) 534 * int m_map 535 * (For internal use only.) 536 * int m_crval 537 * (For internal use only.) 538 * int m_index 539 * (For internal use only.) 540 * int m_indxs 541 * (For internal use only.) 542 * int m_coord 543 * (For internal use only.) 544 * 545 * 546 * Global variable: const char *tab_errmsg[] - Status return messages 547 * ------------------------------------------------------------------ 548 * Error messages to match the status value returned from each function. 549 * 550 *===========================================================================*/ 551 552 #ifndef WCSLIB_TAB 553 #define WCSLIB_TAB 554 555 #ifdef __cplusplus 556 extern "C" { 557 #endif 558 559 560 extern const char *tab_errmsg[]; 561 562 enum tab_errmsg_enum { 563 TABERR_SUCCESS = 0, // Success. 564 TABERR_NULL_POINTER = 1, // Null tabprm pointer passed. 565 TABERR_MEMORY = 2, // Memory allocation failed. 566 TABERR_BAD_PARAMS = 3, // Invalid tabular parameters. 567 TABERR_BAD_X = 4, // One or more of the x coordinates were 568 // invalid. 569 TABERR_BAD_WORLD = 5 // One or more of the world coordinates were 570 // invalid. 571 }; 572 573 struct tabprm { 574 // Initialization flag (see the prologue above). 575 //-------------------------------------------------------------------------- 576 int flag; // Set to zero to force initialization. 577 578 // Parameters to be provided (see the prologue above). 579 //-------------------------------------------------------------------------- 580 int M; // Number of tabular coordinate axes. 581 int *K; // Vector of length M whose elements 582 // (K_1, K_2,... K_M) record the lengths of 583 // the axes of the coordinate array and of 584 // each indexing vector. 585 int *map; // Vector of length M usually such that 586 // map[m-1] == i-1 for coordinate array 587 // axis m and image axis i (see above). 588 double *crval; // Vector of length M containing the index 589 // value for the reference pixel for each 590 // of the tabular coordinate axes. 591 double **index; // Vector of pointers to M indexing vectors 592 // of lengths (K_1, K_2,... K_M). 593 double *coord; // (1+M)-dimensional tabular coordinate 594 // array (see above). 595 596 // Information derived from the parameters supplied. 597 //-------------------------------------------------------------------------- 598 int nc; // Number of coordinate vectors (of length 599 // M) in the coordinate array. 600 int padding; // (Dummy inserted for alignment purposes.) 601 int *sense; // Vector of M flags that indicate whether 602 // the Mth indexing vector is monotonic 603 // increasing, or else decreasing. 604 int *p0; // Vector of M indices. 605 double *delta; // Vector of M increments. 606 double *extrema; // (1+M)-dimensional array of coordinate 607 // extrema. 608 609 // Error handling 610 //-------------------------------------------------------------------------- 611 struct wcserr *err; 612 613 // Private - the remainder are for memory management. 614 //-------------------------------------------------------------------------- 615 int m_flag, m_M, m_N; 616 int set_M; 617 int *m_K, *m_map; 618 double *m_crval, **m_index, **m_indxs, *m_coord; 619 }; 620 621 // Size of the tabprm struct in int units, used by the Fortran wrappers. 622 #define TABLEN (sizeof(struct tabprm)/sizeof(int)) 623 624 625 int tabini(int alloc, int M, const int K[], struct tabprm *tab); 626 627 int tabmem(struct tabprm *tab); 628 629 int tabcpy(int alloc, const struct tabprm *tabsrc, struct tabprm *tabdst); 630 631 int tabcmp(int cmp, double tol, const struct tabprm *tab1, 632 const struct tabprm *tab2, int *equal); 633 634 int tabfree(struct tabprm *tab); 635 636 int tabsize(const struct tabprm *tab, int size[2]); 637 638 int tabprt(const struct tabprm *tab); 639 640 int tabperr(const struct tabprm *tab, const char *prefix); 641 642 int tabset(struct tabprm *tab); 643 644 int tabx2s(struct tabprm *tab, int ncoord, int nelem, const double x[], 645 double world[], int stat[]); 646 647 int tabs2x(struct tabprm *tab, int ncoord, int nelem, const double world[], 648 double x[], int stat[]); 649 650 651 // Deprecated. 652 #define tabini_errmsg tab_errmsg 653 #define tabcpy_errmsg tab_errmsg 654 #define tabfree_errmsg tab_errmsg 655 #define tabprt_errmsg tab_errmsg 656 #define tabset_errmsg tab_errmsg 657 #define tabx2s_errmsg tab_errmsg 658 #define tabs2x_errmsg tab_errmsg 659 660 #ifdef __cplusplus 661 } 662 #endif 663 664 #endif // WCSLIB_TAB 665