[fitswcs] -TAB progress report
Mark Calabretta
mcalabre at atnf.CSIRO.AU
Tue Dec 21 18:36:46 EST 2004
Greetings,
Since WCS Paper III is waiting for the implementation of -TAB, and I
will be on leave till mid-January, I thought a progress report may be
of interest.
-TAB has been implemented, documented, and tested in WCSLIB 3.7
(unreleased) with the exception of the inverse table lookup function,
tabs2x(). The implementation was pretty straightforward and I've only
been able to spend the equivalent of a few full days on it so far. At
this stage I have the following amount of code and documentation:
grus-110% wc tab.h tab.c test/ttab.c
447 2679 19231 tab.h
743 2379 18399 tab.c
123 475 3374 test/ttab.c
1313 5533 41004 total
The WCSLIB header file, tab.h, is appended. Apart from the description
of the arrays in the tabprm struct it looks like most other WCSLIB code.
Comments are invited.
I have yet to start on the FORTRAN wrappers but may defer this till a
later release.
Notes that I made while coding:
* References in Paper III to "indexing array" are somewhat confusing
since these are always one-dimensional, i.e. "indexing vector" would
be a better description.
* In Sect. 6.1.1 it is stated that "The data in the indexing arrays must
be monotonically increasing or decreasing, although two adjacent index
values in the array are allowed to have the same value".
However it should not be considered valid for an index value to appear
more than twice in an index vector, nor for an index value to be
repeated at the start or at the end of an index vector. For example,
the following should be invalid:
2 3 3 3 4 5 ...index 3 repeated more than twice.
8 8 5 3 ...index repeated at start of index vector.
-1 3 3 7 9 9 ...index repeated at end of index vector.
Such repeats in an index vector can never make sense though my
implementation does deal with them except in the case where an index
vector of length > 1 consists solely of a repeated value.
* Eq. (88) is not valid for Psi_k == Psi_{k+1}. The required condition is
either
Psi_k < psi_m <= Psi_{k+1}
or
Psi_k <= psi_m < Psi_{k+1}
for monotonic increasing indexing vectors, and similarly for monotonic
decreasing indexing vectors.
* Paper III gives no hint on how to implement linear interpolation in an
M-dimensional table. I had to resort to the method defined in the
draft of Paper IV.
In testing I was surprised by how inadequate table lookup can be in some
circumstances. However, as I was trying to map coordinates in the
region of a singularity perhaps I expected too much. Nevertheless, it
would be useful to discuss this either in Paper III or IV, probably in
Paper IV.
Mark Calabretta
ATNF
>>>
/*============================================================================
*
* WCSLIB 3.7 - an implementation of the FITS WCS standard.
* Copyright (C) 1995-2004, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library General Public License as published
* by the Free Software Foundation; either version 2 of the License, or (at
* your option) any later version.
*
* This library is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
* Correspondence concerning WCSLIB may be directed to:
* Internet email: mcalabre at atnf.csiro.au
* Postal address: Dr. Mark Calabretta
* Australia Telescope National Facility, CSIRO
* PO Box 76
* Epping NSW 1710
* AUSTRALIA
*
* Author: Mark Calabretta, Australia Telescope National Facility
* http://www.atnf.csiro.au/~mcalabre/index.html
* $Id$
*=============================================================================
*
* WCSLIB 3.7 - C routines that implement tabular coordinate systems as
* defined by the FITS World Coordinate System (WCS) standard. Refer to
*
* "Representations of world coordinates in FITS",
* Greisen, E.W., & Calabretta, M.R. 2002, A&A, 395, 1061 (paper I)
*
* "Representations of spectral coordinates in FITS",
* Greisen, E.W., Valdes, F.G., Calabretta, M.R., & Allen, S.L. 2004, A&A,
* (paper III, in preparation)
*
*
* Summary of routines
* -------------------
* These routines implement the part of the FITS WCS standard that deals with
* tabular coordinates, i.e. coordinates that are defined via a lookup table.
* They define methods to be used for computing tabular world coordinates
* from intermediate world coordinates (a linear transformation of image
* pixel coordinates), and vice versa. They are based on the tabprm struct,
* described in detail below, which contains all information needed for the
* computations. The struct contains some members that must be set by the
* caller, and others that are maintained by these routines, somewhat like a
* C++ class but with no encapsulation.
*
* Three service routines, tabini(), tabcpy(), and tabfree() are provided to
* manage the tabprm struct, and another, tabprt(), to print its contents.
*
* A setup routine, tabset(), computes intermediate values in the tabprm
* struct from parameters in it that were supplied by the caller. The
* struct always needs to be set up by tabset() but it need not be called
* explicitly - see the explanation of tab.flag below.
*
* tabx2s() and tabs2x() implement the WCS tabular coordinate
* transformations.
*
*
* Default constructor for the tabprm struct; tabini()
* ---------------------------------------------------
* tabini() allocates memory for arrays in a tabprm struct and sets all
* members of the struct to default values.
*
* Given:
* alloc int If true, allocate memory for arrays in the tabprm
* struct (see "Memory allocation and deallocation
* below"). Otherwise, it is assumed that pointers to
* these arrays have been set by the caller except if
* they are null pointers in which case memory will be
* allocated for them regardless. (In other words,
* setting alloc true saves having to initalize these
* pointers to zero.)
*
* M int The number of tabular coordinate axes.
* K const int[]
* Vector of length M whose elements (K_1, K_2,... K_M)
* record the lengths of the axes of the coordinate
* array and of each indexing vector. M and K[] are
* used to determine the length of the various tabprm
* arrays and therefore the amount of memory to allocate
* for them.
*
* Given and returned:
* tab struct tabprm*
* Tabular transformation parameters (see below).
* Note that, in order to initialize memory management
* tab->flag should be set to -1 when tab is initialized
* for the first time (memory leaks may result if it had
* already been initialized).
*
* Function return value:
* int Status return value:
* 0: Success.
* 1: Null tabprm pointer passed.
* 2: Memory allocation failed.
* 3: Invalid tabular parameters.
*
*
* Copy routine for the tabprm struct; tabcpy()
* --------------------------------------------
* tabcpy() does a deep copy of one tabprm struct to another, using tabini()
* to allocate memory for its arrays if required. Only the "information to
* be provided" part of the struct is copied; a call to tabset() is required
* to set up the remainder.
*
* Given:
* alloc int If true, allocate memory for the crpix, pc, and cdelt
* arrays in the destination. Otherwise, it is assumed
* that pointers to these arrays have been set by the
* caller except if they are null pointers in which case
* memory will be allocated for them regardless.
* tabsrc const struct tabprm*
* Struct to copy from.
*
* Given and returned:
* tabdst struct tabprm*
* Struct to copy to. tabdst->flag should be set to -1
* if tabdst was not previously initialized (memory
* leaks may result if it was previously initialized).
*
* Function return value:
* int Status return value:
* 0: Success.
* 1: Null tabprm pointer passed.
* 2: Memory allocation failed.
*
*
* Destructor for the tabprm struct; tabfree()
* -------------------------------------------
* tabfree() frees memory allocated for the tabprm arrays by tabini().
* tabini() records the memory it allocates and tabfree() will only attempt
* to free this.
*
* Returned:
* tab struct tabprm*
* Coordinate transformation parameters (see below).
*
* Function return value:
* int Status return value:
* 0: Success.
* 1: Null tabprm pointer passed.
*
*
* Print routine for the tabprm struct; tabprt()
* ---------------------------------------------
* tabprt() prints the contents of a tabprm struct.
*
* Given:
* tab const struct tabprm*
* Tabular transformation parameters (see below).
*
* Function return value:
* int Status return value:
* 0: Success.
* 1: Null tabprm pointer passed.
*
*
* Set up routine; tabset()
* ------------------------
* tabset() allocates memory for work arrays in the tabprm struct and sets up
* the struct according to information supplied within it (see "Tabular
* transformation parameters" below).
*
* Note that this routine need not be called directly; it will be invoked by
* tabx2s() and tabs2x() if the "flag" struct member is anything other than a
* predefined magic value.
*
* Given and returned:
* tab struct tabprm*
* Tabular transformation parameters (see below).
*
* Function return value:
* int Status return value:
* 0: Success.
* 1: Null tabprm pointer passed.
* 3: Invalid tabular parameters.
*
*
* Pixel-to-world transformation; tabx2s()
* ---------------------------------------
* tabx2s() transforms intermediate world coordinates to world coordinates
* using coordinate lookup.
*
* Given and returned:
* tab struct tabprm*
* Tabular transformation parameters (see below).
*
* Given:
* ncoord int The number of coordinates, each of vector length
* nelem int nelem.
* x const double[ncoord][nelem]
* Array of intermediate world coordinates, SI units.
*
* Returned:
* world double[ncoord][nelem]
* Array of world coordinates, in SI units.
* stat int[ncoord]
* Status return value status for each coordinate:
* 0: Success.
* 1: Invalid intermediate world coordinate.
*
* Function return value:
* int Status return value:
* 0: Success.
* 1: Null tabprm pointer passed.
* 3: Invalid tabular parameters.
* 4: One or more of the x coordinates were invalid,
* as indicated by the stat vector.
*
*
* World-to-pixel transformation; tabs2x()
* ---------------------------------------
* tabs2x() transforms world coordinates to intermediate world coordinates.
*
* Given and returned:
* tab struct tabprm*
* Tabular transformation parameters (see below).
*
* Given:
* ncoord int The number of coordinates, each of vector length
* nelem int nelem.
* world const double[ncoord][nelem]
* Array of world coordinates, in SI units.
*
* Returned:
* x double[ncoord][nelem]
* Array of intermediate world coordinates, SI units.
* stat int[ncoord]
* Status return value status for each vector element:
* 0: Success.
* 1: Invalid world coordinate.
*
* Function return value:
* int Status return value:
* 0: Success.
* 1: Null tabprm pointer passed.
* 3: Invalid tabular parameters.
* 5: One or more of the world coordinates were
* invalid, as indicated by the stat vector.
*
*
* Tabular transformation parameters
* ----------------------------------
* The tabprm struct consists of the following elements that must be
* supplied:
*
* int flag
* This flag must be set to zero whenever any of the following tabprm
* structure members are set or changed. This signals the
* initialization routine, tabset(), to recompute intermediaries.
*
* int M
* Number of tabular coordinate axes.
*
* int *K
* Pointer to the first element of a vector of length M whose elements
* (K_1, K_2,... K_M) record the lengths of the axes of the coordinate
* array and of each indexing vector.
*
* int *map
* Pointer to the first element of a vector of length M that defines
* the association between axis m in the M-dimensional coordinate array
* (1 <= m <= M) and the indices of the intermediate world coordinate
* and world coordinate arrays, x[] and world[], in the argument lists
* for tabx2s() and tabs2x().
*
* When x[] and world[] contain the full complement of coordinate
* elements in image-order, as will usually be the case, then
* map[m-1] == i-1 for axis i in the N-dimensional image (1 <= i <= N).
* In terms of the FITS keywords
*
* map[PVi_3a - 1] == i - 1.
*
* However, a different association may result if x[], for example,
* only contains a (relevant) subset of intermediate world coordinate
* elements. For example, if M == 1 for an image with N > 1, it is
* possible to fill x[] with the relevant coordinate element with nelem
* set to 1. In this case map[0] = 0 regardless of the value of i.
*
* double *crval
* Pointer to the first element of a vector of length M whose elements
* contain the index value for the reference pixel for each of the
* tabular coordinate axes.
*
* double *coord
* Pointer to the first element of the tabular coordinate array,
* treated as though it were defined as
*
* double coord[K_M]...[K_2][K_1][M]
*
* i.e. with the M dimension varying fastest so that the M elements of
* a coordinate vector are stored contiguously in memory.
*
* double **index
* Pointer to the first element of a vector of length M of pointers to
* indexing vectors of lengths (K_1, K_2,... K_M).
*
* The remaining members of the tabprm struct are maintained by tabset() and
* must not be modified elsewhere:
*
* int *sense
* Pointer to the first element of a vector of length M whose elements
* indicate whether the corresponding indexing vector is monotonic
* increasing (+1), or decreasing (-1).
*
* double *p0, *delta
* Pointer to the first element of a vector of length M of interpolated
* indices into the coordinate array such that Upsilon[m], as defined
* in Paper III, is equal to p0[m] + delta[m].
*
*
* Vector arguments
* ----------------
* Arrays of intermediate world coordinates and world coordinates are two
* dimensional, i.e. x[ncoord][nelem].
*
* Note that the function prototypes must declare two-dimensional arrays as
* one-dimensional to avoid compiler warnings about declaration of
* "incomplete types". This was considered preferable to declaring them as
* simple pointers-to-double which gives no indication that storage is
* associated with them.
*
*
* Memory allocation and deallocation
* ----------------------------------
* tabini() allocates memory for the K, map, crval, coord, and index arrays
* (including the arrays referenced by index[]) in the tabprm struct. It is
* provided as a service routine; usage is optional, and the caller is at
* liberty to set these pointers independently.
*
* tabset() also allocates memory for the sense, p0, and delta arrays. The
* caller must not modify these.
*
* tabini() maintains a record of memory it has allocated and this is used
* by tabfree() which tabini() uses to free any memory that it may have
* allocated on a previous invokation. Thus it is not necessary for the
* caller to invoke tabfree() separately if tabini() is invoked repeatedly on
* the same tabprm struct. Likewise, tabset() deallocates memory that it
* may have allocated in the same tabprm struct on a previous invokation.
*
* However, a memory leak will result if a tabprm struct goes out of scope
* before the memory has been free'd, either by tabfree() or otherwise.
* Likewise, if the tabprm struct itself has been malloc'd and the allocated
* memory is not free'd when the memory for the struct is free'd. A leak may
* also arise if the caller interferes with the array pointers in the
* "private" part of the tabprm struct.
*
* Beware of making a shallow copy of a tabprm struct by assignment; any
* changes made to allocated memory in one would be reflected in the other,
* and if the memory allocated for one was free'd the other would reference
* unallocated memory. Use tabcpy() instead to make a deep copy.
*
*
* Status return values
* --------------------
* Error messages to match the status value returned from each function are
* are encoded in the tab_errmsg character array.
*
*
* Accuracy
* --------
* No warranty is given for the accuracy of these routines (refer to the
* copyright notice above); intending users must satisfy for themselves their
* adequacy for the intended purpose. However, closure effectively to within
* double precision rounding error was demonstrated by test routine ttab.c
* which accompanies this software.
*
*===========================================================================*/
#ifndef WCSLIB_TAB
#define WCSLIB_TAB
#ifdef __cplusplus
extern "C" {
#endif
extern const char *tab_errmsg[];
#define tabini_errmsg tab_errmsg
#define tabcpy_errmsg tab_errmsg
#define tabfree_errmsg tab_errmsg
#define tabprt_errmsg tab_errmsg
#define tabset_errmsg tab_errmsg
#define tabx2s_errmsg tab_errmsg
#define tabs2x_errmsg tab_errmsg
struct tabprm {
/* Initialization flag (see the prologue above). */
/*-----------------------------------------------------------------------*/
int flag; /* Set to zero to force initialization. */
/* Parameters to be provided (see the prologue above). */
/*-----------------------------------------------------------------------*/
int M; /* Number of tabular coordinate axes. */
int *K; /* Vector of length M whose elements */
/* (K_1, K_2,... K_M) record the lengths of */
/* the axes of the coordinate array and of */
/* each indexing vector. */
int *map; /* Vector of length M usually such that */
/* map[m-1] == i-1 for coordinate array */
/* axis m and image axis i (see above). */
double *crval; /* Vector of length M containing the index */
/* value for the reference pixel for each */
/* of the tabular coordinate axes. */
double *coord; /* (1+M)-dimensional tabular coordinate */
/* array (see above). */
double **index; /* Vector of pointers to M indexing vectors */
/* of lengths (K_1, K_2,... K_M). */
/* Information derived from the parameters supplied. */
/*-----------------------------------------------------------------------*/
int *sense; /* Vector of M flags that indicate whether */
/* the Mth indexing vector is monotonic */
/* increasing, or else decreasing. */
int *p0; /* Vector of M indices. */
double *delta; /* Vector of M increments. */
int m_flag, m_M, m_N; /* The remainder are for memory management. */
int *m_K, *m_map;
double *m_crval, *m_coord, **m_index;
int set_M;
};
#define TABLEN (sizeof(struct tabprm)/sizeof(int))
int tabini(int, int, const int[], struct tabprm *);
int tabcpy(int, const struct tabprm *, struct tabprm *);
int tabfree(struct tabprm *);
int tabprt(const struct tabprm *);
int tabset(struct tabprm *);
int tabx2s(struct tabprm *, int, int, const double[], double[], int[]);
int tabs2x(struct tabprm *, int, int, const double[], double[], int[]);
#ifdef __cplusplus
};
#endif
#endif /* WCSLIB_TAB */
More information about the fitswcs
mailing list