[fitswcs] -TAB progress report

[Mark Calabretta]

On Wed 2004/12/22 10:36:46 +1100, Mark Calabretta wrote
in a message to: fitswcs at nrao.edu

>-TAB has been implemented, documented, and tested in WCSLIB 3.7
>(unreleased) with the exception of the inverse table lookup function,

Greetings,

-TAB is now completely implemented in WCSLIB 4.0.  Here are a few notes 
I made in the process:

* Inversion of 1-dimensional tables is implemented via a separate
  algorithm based on Eq. (89) that is much simpler and faster than that
  for multi-dimensional tables.

* Repeated index values slightly complicate the inversion of a 1-D
  table.  This is not apparent in the prescription given by Eq. (89).

  For the header in Table 11, the time computed for pixel coordinate
  p4 = 2.0 (psi_4 = 1.5) is 1993.28454, being halfway between 1993.28451
  and 1993.28456.

  To invert that using Eq. (89) you would scan the coordinate vector
  from the start looking for a pair of values that encompass 1993.28454,
  namely 1997.84631 and 1993.28451.  Applying Eq. (89) gives pixel
  coordinate 1.5 (psi_4 = 1.0) - in disagreement with the original
  value.

  In order to resolve this non-closure you need to look at the index
  vector and check that the two indices (Psi_k) differ for the two
  coordinate values, 1997.84631 and 1993.28451.  In this case they are
  equal so the solution must be rejected.  The correct solution is
  obtained from the next pair of coordinates.

* Implementation of the general multi-dimensional inverse table lookup
  is more difficult than supposed.  The table defines a function of M
  variables whose value is a vector of length M, where M is not known in
  advance; inverting this is much more complicated than inverting a
  scalar-valued function.  I have implemented a compact interval-halving
  algorithm based on recursion but it's difficult to see how to do it in
  a reasonable way without recursion.  

* For the table used in WCSLIB test program ttab3 (size 235 x 271 x 2 =~
  125000) the inverse table lookup is rather slow because of the table
  searching required.  Caching the maxima and minima for each row can
  speed it up, though the speedup is very dependent on the exact nature
  of the coordinate system described by the table.

* As currently defined, the pixel-to-world calculation is fast, while
  world-to-pixel is slow.  This is probably the right way round; pixel-
  to-world is used for reporting cursor coordinates, while world-to-
  pixel is used for drawing graticules.  However, slow graticules may
  present a problem for zooming and panning.

* Blanking (undefined coordinate values) should be considered; the
  interior of a voxel that has a blank in one corner would be considered
  to be blank.  This could potentially speed up the inverse table
  lookup.

* Repeated index values make little sense for M > 1 and should not be
  allowed.  They would make the implementation for M > 1 too difficult.

* Km == 1 (i.e. a degenerate table axis) causes all sorts of problems in
  coding, perhaps it should be disallowed.

WCSLIB test program ttab1 verifies closure for a 1-D and a 2-D table.

ttab2 demonstrates the nature of linear interpolation in 2-D; the level
curves are either parallel, equi-spaced straight lines (i.e. a plane
surface), or else hyperbolic arcs.

ttab3 constructs a 2-D table that approximates Bonne's projection and
then uses inverse table lookup to draw a graticule.  The result is a
very good approximation of the true graticule.

Mark Calabretta
ATNF