From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001
From: Joseph Hunkeler <jhunkeler@gmail.com>
Date: Wed, 8 Jul 2015 20:46:52 -0400
Subject: Initial commit

---
 sys/mwcs/MWCS.hlp | 1026 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 1026 insertions(+)
 create mode 100644 sys/mwcs/MWCS.hlp

(limited to 'sys/mwcs/MWCS.hlp')

diff --git a/sys/mwcs/MWCS.hlp b/sys/mwcs/MWCS.hlp
new file mode 100644
index 00000000..4f77144f
--- /dev/null
+++ b/sys/mwcs/MWCS.hlp
@@ -0,0 +1,1026 @@
+.help MWCS Oct89 "Mini-WCS Interface"
+
+.ce
+\fBMini-WCS Interface\fR
+.ce
+Doug Tody
+.ce
+October 1989
+
+
+.nh
+Introduction
+
+    The mini-WCS interface represents a first cut at the general problem
+of representing a linear or nonlinear world coordinate system (WCS).
+While some of the harder problems are avoided and the general WCS problem
+remains to be solved, the current interface should be largely upwards
+compatible with future versions of the interface.  The main items omitted
+from this initial version of the interface are support for general nonlinear
+world coordinate systems, particularly support for modeling of geometric
+distortions and arbitrary application defined coordinate mapping functions.
+Limited support is provided for the projective geometries.
+
+.nh
+WCS Definition
+.nh 2
+Linear Transformations
+
+    Any linear transformation consisting of some combination of a shift,
+rotate, axis-flip, scale change, etc. can be expressed as
+
+.ks
+.nf
+	|x'|   |a b|   |x|   |u|
+	|  | = |   | * | | + | |                                   [2.1]
+	|y'|   |c d|   |y|   |v|
+.fi
+.ke
+
+where [x,y] are the input coordinates, [x',y'] are the transformed
+coordinates, [a,b,c,d] is a rotation matrix, and [u,v] is a shift vector.
+
+For example, the X term of a combination of a rotation about a point [x0,y0]
+plus a shift to an offset [x1,y1] may be expressed as
+
+.ks
+.nf
+	x' = a(x - x0) + b(y - y0) + x1
+	   = ax - ax0 + by - by0 + x1
+	   = ax + by + u
+.fi
+.ke
+whence
+.nf
+	 u = x1 - ax0 - by0
+and                                                                [2.2]
+	 v = y1 - cx0 - dy0
+.fi
+
+Another way of expressing this is to note that [U,V] is the transform
+of the origin [x,y]=[0,0] of the original coordinate system.
+There is nothing special about the rotation point; a rotation about
+any point [x,y] is equivalent to a rotation about the origin followed
+by a translation equal to the distance of the rotation point from the origin.
+
+The inverse transformation is given by
+
+.nf
+	|x|   |  -1|    / |x'|   |u| \
+	| | = | A  | * <  |  | - | |  >                            [2.3]
+	|y|   |    |    \ |y'|   |v| /
+.fi
+
+where A**(-1) is the inverse of the rotation matrix [a,b,c,d].
+
+.nh 2
+World Coordinate Systems
+
+    A world coordinate system (WCS) defines the transformation between
+a physical coordinate system (e.g., pixel coordinates in a reference image),
+and world coordinates expressed in some arbitrary units.  A two dimensional
+WCS can be expressed as
+
+.ks
+.nf
+	(x',y') = F (l,m, Wx,Wy)
+	(l,m) = [CD] * (x-Rx, y-Ry)                                [2.4]
+
+where
+
+	x,y	Are the coordinates of a point in the physical system.
+	l,m	Is a linearly transformed representation of the point.
+	x',y'	Are the coordinates of the same point in the world system.
+	F	Is the WCS function, possibly a nonlinear function.
+	Rx,Ry	Define the reference point in the physical system.
+	Wx,Wy	Are the world coordinates of the reference point.
+	[CD]	Is the coefficient determination (CD) matrix.
+.fi
+.ke
+
+The notation [CD]*(x,y) denotes a matrix multiply of the CD matrix [CD] and
+the vector (x,y), i.e., a linear transformation of the vector (x,y).
+If the WCS contains a nonlinear component, as for a sky projection,
+this is described by the function F in terms of the intermediate
+coordinates (l,m), e.g., the displacement in degrees from the reference
+point.  Separation of the WCS into linear and nonlinear components allows
+full specification of linear systems using only the basic interface,
+and simplifies the representation of the nonlinear part of the WCS.
+The nonlinear component itself (F) may be an object of arbitrary complexity.
+
+In the case of a simple 2D linear WCS with no rotation this reduces to
+
+.ks
+.nf
+	x' = (x - Rx) * CD[1,1] + Wx
+	y' = (y - Ry) * CD[2,2] + Wy
+.fi
+.ke
+
+In the general case the world system may be rotated with respect to the
+physical system (original image matrix), hence the WCS must include a
+rotation term.  Specifying this as a general linear transformation expressed
+as a matrix multiplication allows the representation of such transformations
+as conversion between skewed and cartesian coordinates as well as the
+more conventional rotation and scale transformation.
+
+The CD matrix representation (developed by STScI and now also associated
+with FITS), in addition to allowing specification of a linear transformation,
+is responsible for converting between the coordinate units used in the
+physical system and those used in the world system (more precisely,
+in the general case the units of (l,m) may differ from those of the
+world system, since F can also change the units).
+It is even possible for the world system to use different units on
+different axes, so long as a rotation is not defined between axes with
+different units.
+
+For example, if the WCS is used to describe an image cube, following
+application of the CD matrix axes 1 and 2 might have units of arc seconds,
+and axis 3 frequency.  In this case rotation would be defined only between
+axes 1 and 2, i.e., the off-diagonal CD matrix terms CD[1:2,3] and CD[3,1:2]
+must be zero, with CD[3,3] giving the scale for axis 3 independently of any
+rotation between axes 1 and 2.  This restriction on rotation between
+dissimilar axes applies only to the world system described by the CD matrix.
+As we shall see in the next section, when the WCS refers to an image,
+arbitrary rotations of the raw pixel matrix are still possible by using
+a separate pixel space transformation to describe transformations
+of the image matrix.
+
+.nh 3
+WCS Rotation Between Dissimilar Axes
+
+    To see why rotation between dissimilar axes is disallowed in some
+circumstances, note that the CD matrix, since it combines a rotation
+(or other pixel space linear transformation) and units conversion in one
+operation, can be expressed as follows in the case of a two dimensional system.
+
+The CD matrix is used as follows:
+
+.ks
+.nf
+	| l |   |    |   | x |
+	|   | = | CD | * |   |
+	| m |   |    |   | y |
+.fi
+.ke
+
+The CD matrix is constructed as follows:
+
+.ks
+.nf
+	|    |   | Dx   0 |   | a  b |
+	| CD | = |        | * |      |
+	|    |   | 0   Dy |   | c  d |
+.fi
+.ke
+
+where (Dx,Dy) is the units conversion matrix, and (a,b,c,d) is the
+rotation matrix.  This is a completely general representation, i.e.,
+any linear transformation may be specified by the matrix (a,b,c,d)
+and combined with the units conversion matrix, since the rotation
+matrix rotates the physical system to align the axes with those of
+the world coordinate system.
+
+The problem comes if we try to \fIrotate the CD matrix\fR.
+Although the CD matrix can express any
+rotation between the physical and world system, once the CD matrix
+has been formed the units conversion and rotation matrices cannot
+be recovered.  In general, further rotation of the system described by
+the CD matrix requires that we rotate the matrix (a,b,c,d), rather than
+the CD matrix itself.  The only exception occurs when Dx=Dy (similar axes),
+in which case the CD matrix and rotation matrix are equivalent except
+for a constant.  Hence, rotations between dissimilar axes of the system
+described by the (already formed) CD matrix are disallowed.  A special
+case is rotation is some multiple of 90 degrees, which can be represented
+by an axis swap or flip.
+
+.nh 2
+MWCS Coordinate System Representation
+
+    The coordinate system representation used by the MWCS interface consists of
+two components called the \fBLterm\fR and \fBWterm\fR, specifying independent
+logical and world transformations relative to a physical, cartesian coordinate
+system.  Three types of coordinate systems are defined, as outlined below.
+
+.ls
+.ls PHYSICAL
+The physical coordinate system is the raw coordinate system of the data.
+In the case of an image, the physical coordinate system refers to the pixel
+coordinates of the original data frame.  All other coordinates systems are
+defined in terms of the physical system (reference frame).
+.le
+.ls LOGICAL
+The logical coordinate system is defined by the \fILterm\fR in terms of the
+physical coordinate system.  In the case of an image, the logical coordinate
+system specifies raw pixel coordinates relative to some image section or
+derived image, i.e., the coordinates used for image i/o.  In the MWCS the
+Lterm specifies a simple linear transformation, in pixel units, between
+the original physical image matrix and the current image section.
+.le
+.ls WORLD
+The world coordinate system is defined by the \fIWterm\fR in terms of the
+physical coordinate system.  Any number of different kinds of world coordinate
+systems are conceivable.  Examples are the tangent (gnonomic) projection,
+specifying right ascension and declination relative to the original data
+image, or any linear WCS, e.g., a linear dispersion relation for spectral
+data.  Multiple world coordinate systems may be simultaneously defined in
+terms of the same physical system.
+.le
+.le
+
+The following observations apply to the behavior of MWCS as applied
+to image data.
+.ls
+.ls 4 1.
+Any linear transformation of the image matrix (shift, scale change,
+axis flip, etc.) affects only the Lterm.  The revised MWCS for the new
+image or image section may be computed merely by doing a linear
+transformation of the Lterm.
+.le
+.ls 4 2.
+If multiple world coordinate systems are associated with an image,
+all share the same Lterm.
+.le
+.ls 4 3.
+Geometric distortion of an image (not currently supported by MWCS) is
+a pixel space operation, i.e. a generalization of the Lterm, hence is
+independent of the WCS.
+.le
+.le
+
+In general, the physical and world coordinate systems are defined whenever
+a new image is created, e.g., by a task such as RFITS.  A new-copy type
+operation, such as most transformations performed by IMAGES tasks, affects
+only the Lterm.
+
+Although we normally speak in terms of images, MWCS is not
+limited to applications involving images.  For example, the physical
+coordinate system could just as well be a graphics frame buffer, and the
+logical coordinate system a pixrect.  A greyscale transformation is
+an example of a non-image WCS.  MWCS, or the successor interface,
+will eventually be used in GIO, e.g., for cursor readback.
+
+Since the Wterm includes the CD matrix, which defines a linear
+transformation, and linear transformations can be combined, in principle
+it is possible to combine the Lterm and Wterm to define a single
+transformation from logical to world coordinates.  In practice this
+can run into problems, as not all pixel space rotations may be
+representable by the CD matrix (since the latter may define different
+world space units on different axes).  Furthermore, if multiple WCS
+are defined, and the WCS are defined in terms of the logical system,
+it would be inefficient to have to transform each WCS to the new logical
+system each time a linear transformation of the data is performed
+(e.g., every time an image is opened with an image section).
+
+For images, the common coordinate transformations are image section (logical)
+coordinates to world coordinates and vice versa, and section coordinates
+to physical coordinates and vice versa.  The physical coordinate system
+can be regarded as a special case of a world coordinate system (the "pixel"
+coordinate system) defined relative to logical image section coordinates.
+
+For example, to convert IMIO image coordinates to world coordinates,
+the interface will first apply the inverse of the Lterm to determine the
+coordinates in the physical system, then apply the Wterm for the desired WCS
+to compute the world coordinates.  If the Wterm is linear and the Lterm
+does not define any rotations between dissimilar axes, the two operations
+can be combined for a more efficient coordinate transformation.
+
+An arbitrary number of world coordinate systems may be defined over the
+same domain in the physical coordinate system.  Every WCS has a name
+uniquely specifying the WCS type.  A WCS may also have attributes such as
+units, axis labels, and numeric output formats specified independently for
+each axis, as well as arbitrary user defined WCS attributes.
+
+.nh 3
+Lterm Representation
+
+    The Lterm is defined by the terms of a general linear transformation,
+as shown in equation [2.1].  For example, in the case of a 2D system the
+following quantities must be given to define the Lterm.
+
+.ks
+.nf
+	[CD] = [a,b,c,d]	rotation matrix
+	tv[] = [u,v]		translation vector
+.fi
+.ke
+
+This defines the transformation between the physical and logical coordinate
+systems, i.e., applying the transformation to a pair of physical coordinates
+[x,y] yields the corresponding logical coordinates [x',y'].  MWCS will
+automatically compute the inverse transformation when asked to convert
+between logical coordinates and physical or world coordinates.
+
+.nh 3
+Wterm Representation
+
+    The Wterm defines the transformation between the physical system and
+some arbitrary world coordinate system.  The Wterm is defined by the
+following quantities:
+
+.nf
+	R[]		reference coordinates in physical system
+	W[]		world coordinates at the reference point
+	[CD]		coordinate determination matrix
+	wtype		type of WCS (function name string)
+	wattr		WCS attributes (string, opaque outside interface)
+.fi
+
+The point R, also known as the \fIreference pixel\fR when dealing
+with image data, defines the origin of the world coordinate system in
+the physical system.  The world coordinates at the reference point are
+given by the vector W (at least for a linear WCS; in general the meaning
+of the W term depends upon the WCS type).  The CD matrix defines any
+rotation between the physical and world systems, as well as the scale
+conversion needed to convert between physical and world (or linear world)
+coordinates.
+
+Although the function name or type \fIwtype\fR is
+accessible to applications, the details of what the WCS means, and how
+it is evaluated, are intended to be internal to the interface, hence
+the use of strings to pass in the WCS information.  Functions complex
+enough to require coefficients should pass the extra information in via
+the \fIwattr\fR term.  WCS attributes such as the axis units and labels
+are also passed in via \fIwattr\fR.
+
+In the general case there may be any number of different WCS types.
+In the case of MWCS, however, only a predefined set of WCS types are
+supported, since the code for each WCS is wired into the interface.
+The predefined WCS types (as selected by \fIwtype\fR) are the following.
+
+.ks
+.nf
+	\fIWtype\fR			\fIDescription\fR
+
+	linear			simple linear WCS
+	sampled			sampled WCS function
+	(TAN etc.)		the sky projections
+.fi
+.ke
+
+A \fIlinear\fR WCS is specified by the physical and world coordinates of the
+reference point, and the row or rows of the CD matrix pertaining to the axes
+to which the WCS is assigned.  A linear WCS is completely specified by the
+linear term of the standard WCS representation.
+  
+A \fIsampled\fR WCS is specified by an array of (physical, world) coordinate
+pairs, i.e., an array of reference points, sampling the linear WCS function
+for that axis.  In the limiting case, for sampled (pixel) data, there is one
+(physical, world) point on the WCS curve for each data point.  If the WCS
+function is smooth a coarser sampling can be used to approximate the curve,
+using some form of interpolation to evaluate the function.  In the MWCS,
+a sampled function must be one-dimensional, i.e., associated with a single
+axis (higher dimensional surfaces can be represented so long as the axes
+are independent).
+
+Note that the sampled function is expressed in terms of the \fIoffset\fR
+from the reference point in both the physical and world systems.
+Any analytic function, e.g., polynomial or spline, can be sampled and
+later reconstructed from the sampled curve with no significant error,
+provided the function type and order are known and there are sufficient
+sample points to determine the system.
+The advantage of the sampled representation is that it is independent of
+the function type, and can be used to fit any analytic function when the
+time comes to evaluate the curve.
+
+The sky projections are a special case used with astronomical direct images. 
+The principal example is the gnomonic projection, the projection of the
+celestial sphere onto a plane tangent at the reference point.
+The sky projections are completely specified by the reference point and the
+standard CD matrix, plus the WCS name which specifies the type of projection,
+e.g., "gnomonic", "sine", "arc", and so on.
+
+The WCS attributes which can be set by the \fIwattr\fR string consist of
+a number of standard attributes, plus an arbitrary number of additional
+WCS specific attributes.  Examples of standard attributes include "system",
+"wtype", "units", "label", etc.  A list of the standard WCS attributes is
+given in section 3.3.1.
+
+.nh
+Interface Overview
+
+    The MWCS interface is a stand-alone interface implementing the linear
+and world coordinate transformation abstractions.  While the interface
+is designed with the typical application to image data in mind, MWCS is
+intended as a general coordinate transformation facility for use with any
+type of data, as an embedded interface in other software, including system
+interfaces such as IMIO and GIO as well as user applications.
+
+.nh 2
+Object Creation and Storage
+
+    The MWCS interface routines used to create or access MWCS objects, or
+save and restore MWCS objects in external storage, are summarized below.
+
+.nf
+		   mw = mw_open (bufptr|NULL, ndim)
+		 mw = mw_openim (im)
+		mw = mw_newcopy (mw)
+		       mw_close (mw)
+
+			mw_load (mw, bufptr)
+		  len = mw_save (mw, bufptr, buflen)
+	       mw_[load|save]im (mw, im)
+.fi
+
+A new MWCS object, initialized either to a unitary transformation of
+dimension \fIndim\fR or to the encoded MWCS in the input buffer,
+is created with \fImw_open\fR.  A MWCS object is be
+created and initialized from an image with \fImw_openim\fR; if the referenced
+image does not currently have any WCS information associated with it,
+a unitary pixel WCS will be created.  The \fImw_newcopy\fR operation
+creates a new MWCS object as a copy of an existing one, as one might wish
+to do prior to modifying a WCS.  When a descriptor is no longer needed it
+should be returned with \fImw_close\fR.
+
+A MWCS object (descriptor) is a memory object.  To encode a MWCS in an opaque
+machine independent binary array, e.g., for storage in a file or transmission
+through a datastream, \fImw_save\fR is called with the \fIchar\fR pointer of
+the buffer in which the encoded MWCS is to be placed.  If the buffer pointer is
+NULL a buffer will be created and the pointer returned, and if a valid buffer
+is passed it will be resized as necessary to store the encoded object.
+An encoded MWCS object is reloaded into a descriptor with \fImw_load\fR.
+A MWCS may be stored or updated in an image header with \fImw_saveim\fR,
+or loaded from the image header into a descriptor with \fImw_loadim\fR.
+These are the only interface routines with knowledge of the parameter names,
+etc., used to store WCS information in image headers.
+
+.nh 2
+Coordinate Transformation Procedures
+
+    The MWCS procedures used to perform coordinate transformations,
+and to modify or examine the Lterm and Wterm, are summarized below.
+
+.nf
+		 ct = mw_sctran (mw, system1, system2, axes)
+	  ndim = mw_gctran[r|d] (ct, ltm, ltv, axtype1, axtype2, maxdim)
+		      mw_ctfree (ct)
+
+	    x2 = mw_c1tran[r|d] (ct, x1)
+		 mw_v1tran[r|d] (ct, x1, x2, npts)
+		 mw_c2tran[r|d] (ct, x1,y1, x2,y2)
+		 mw_v2tran[r|d] (ct, x1,y1, x2,y2, npts)
+		  mw_ctran[r|d] (ct, p1, p2, ndim)
+		  mw_vtran[r|d] (ct, v1, v2, ndim, npts)
+.fi
+
+The procedures \fImw_[cv][12]tran[rd]\fR perform coordinate
+transformations for individual coordinates or coordinate vectors,
+for one or two dimensional systems, for coordinates of type real or double.
+The general N dimensional case is handled by the \fImw_[cv]tran[rd]\fR
+procedures, which transform \fIndim\fR-dimensional points (\fImw_ctran\fR)
+or point vectors (\fImw_vtran\fR).  A single point is specified as a
+vector of length \fIndim\fR; a point vector is expressed as an array of
+points, i.e., a 2-dimensional array V[I,J], where the index I refers to
+the axis within a point vector, and where the index J refers to the point.
+The notation V[1,j] references the 1-dimensional point vector for point J.
+
+The direction of the transformation, and the axes for which the transformation
+is to be performed, is determined by a prior call to \fImw_sctran\fR,
+which specifies the input and output coordinate systems, and performs the
+initialization necessary for efficient evaluation of a series of
+transformations.  A pointer to the optimized transformation descriptor is
+returned, to allow two or more transformations to be prepared and used
+simultaneously without having to repeat the setup, which can be considerably
+more expensive than coordinate evaluation.  The transformation descriptor
+should be freed when no longer needed, else it will be freed automatically
+when the MWCS is closed.
+
+A coordinate system is specified to \fImw_sctran\fR by its name.
+The following standard systems are predefined.
+Additional WCS names may be defined by the application.
+
+.ks
+.nf
+	"logical"		The logical system
+	"physical"		The physical system
+	"world"			The default world system
+	(user-wcs)		User defined systems
+.fi
+.ke
+
+Strings are used to specify the coordinate systems in order to allow user
+defined and named systems to be added at runtime.
+The use of a setup procedure to specify the desired transformation allows
+new types of coordinate transformations to be easily added, for example
+mixed conversions, as for a 2-dimensional system where the X and Y components
+of a coordinate pair belong to different coordinate systems, or computation
+of the derivative at a point.  In MWCS, only simple conversions between any
+two of the physical, logical, and world coordinate systems are supported.
+
+Specification of the axes for which the coordinate transformation is desired
+is necessary for the more complex systems, since there may be different,
+often quite independent coordinate systems defined on different axes.
+The axes for which the transformation is to be prepared are specified as
+a bitmask.  The default, if the mask is zero, is to use axes starting with
+1, up to the number required to satisfy the given dimension transformation.
+
+For example, to convert two dimensional image coordinates (section relative)
+to world coordinates in the default WCS:
+
+.nf
+	call mw_sctran (mw, "logical", "world", 3B)
+	call mw_c2tranr (mw, px,py, wx,wy)
+.fi
+
+Multiple independent world coordinate systems may be defined relative to
+the same physical system.  Most applications, however, are best written
+as if there were only one world system, with the coordinate system to be
+used being switched about transparently to the application.  For this
+reason there is no WCS number argument to the MWCS procedures, and
+the "world" system specifies the \fIcurrent default\fR WCS.  If a MWCS object
+defines multiple world coordinate systems, a \fImw_ssystem\fR call is used
+to select the WCS to be used.  This could be used, for example, to change
+the units appearing on plots in a graphics application, transparently
+to the application.
+
+.nh 2
+Coordinate System Specification
+
+    The MWCS procedures used to enter, modify, or inspect the MWCS
+logical and world coordinate transformations are summarized in the figure
+below.
+
+The procedures \fImw_[sg]lterm\fR are used to directly enter
+or inspect the Lterm, which consists of the linear transformation matrix
+\fIltm\fR and the translation vector \fItv\fR, both of dimension \fIndim\fR,
+defining the transformation from the physical system to the logical system.
+If the logical system undergoes successive linear transformations,
+\fImw_translate\fR may be used to translate, rather than replace,
+the current Lterm, where the given transformation matrix and translation
+vector refer to the relative transformation undergone by the logical system.
+This will always work since the Lterm is initialized to the identity matrix
+when a new MWCS object is created.  The routines \fImw_rotate\fR,
+\fImw_scale\fR, and \fImw_shift\fR provide a convenient front-end to
+\fImw_translate\fR for the more common types of translations.
+
+Specification of the Wterm is somewhat more complicated.  The Wterm for
+a new WCS should first be created and initialized for a system of the given
+dimensionality with \fImw_newsystem\fR.  The linear portion of the Wterm,
+i.e., the CD matrix and the coordinates of the reference point
+in the physical and world systems, and the WCS dimension, may then be entered
+with \fImw_swterm\fR and queried with \fImw_gwterm\fR.
+
+.nf
+	     mw_[s|g]lterm[r|d] (mw, ltm, ltv, ndim)
+	      mw_translate[r|d] (mw, ltv_1, ltm, ltv_2, ndim)
+		      mw_rotate (mw, theta, center, axes)
+		       mw_scale (mw, scale, axes)
+		       mw_shift (mw, shift, axes)
+
+		   mw_newsystem (mw, system, ndim)
+		 mw_[s|g]system (mw, system[, maxch])
+		  mw_[s|g]axmap (mw, axno, axval, ndim)
+		    mw_bindphys (mw)
+
+	     mw_[s|g]wterm[r|d] (mw, r, w, cd, ndim)
+		      mw_swtype (mw, axis, naxes, wtype, wattr)
+	     mw_[s|g]wsamp[r|d] (mw, axis, pv, wv, npts)
+		 mw_[s|g]wattrs (mw, axis, attribute, valstr[, maxch])
+.fi
+
+The world portion of the Wterm is unusual in that the type of WCS may be
+specified independently for each axis.  The WCS function type \fIwtype\fR,
+and any attributes \fIwattr\fR, are specified for the indicated \fIaxes\fR
+with \fImw_swtype\fR.  The axes specified are those required to evaluate
+the named function.
+
+In the case of an axis of type "sampled", the sampled WCS function must
+also be entered via a call to \fImw_swsamp\fR, and may later be retrieved
+with \fImw_gwsamp\fR.  The WCS function is defined as an \fIoffset\fR from
+the reference point in both the physical and world systems, e.g.,
+the vector \fIWv\fR will be added to the world coordinates
+produced by interpolating the sampled function (this can of course be
+defeated by setting R or W to zero).
+
+A WCS always has a number of predefined \fIattributes\fR, and may also
+have any number of user defined, or WCS specific, attributes.  These are
+defined when the WCS is created, in the \fIwattr\fR argument input to
+\fImw_swtype\fR, or in a subsequent call to \fImw_swattrs\fR.  The WCS
+attributes for a specific axis may be queried with the function
+\fImw_gwattrs\fR.  Attribute values may be modified, or new attributes defined,
+with \fImw_swattrs\fR.  The issue of WCS attributes is discussed further
+in the next section.
+
+.nh 3
+WCS Types and Attributes
+
+    The WCS attributes which can be set by the \fIwattr\fR term consist of
+a number of standard attributes, plus an arbitrary number of additional
+WCS specific (application defined) attributes.  The following standard
+attributes are reserved (but not necessarily defined) for each WCS:
+
+.nf
+	"units"		axis units ("pixels", etc.)
+	"label"		axis label, for plots
+	"format"	axis numeric format, for tick labels
+	"wtype"		WCS type, e.g., "linear"
+.fi
+
+In addition, the following are defined for the entire WCS,
+regardless of the axis:
+
+.nf
+	"system"	system name (logical, physical, etc.)
+	"object"	external object with which WCS is associated
+.fi
+
+For example, to determine the WCS type for axis 1:
+
+	call mw_gwattrs (mw, 1, "wtype", wtype, SZ_WTYPE)
+
+The (world coordinate) system name \fIsystem\fR is what is used, e.g.,
+to select a WCS in a call to \fImw_ssystem\fR, or define a coordinate
+transformation in a call to \fImw_sctran\fR.  Note that the system name
+"world" is actually only an alias for the \fIdefault world system\fR.
+This may be any primary system, i.e., the logical or physical system,
+or a user defined world system.  The initial default world system may
+be specified by the user by predefining the environment variable
+\fBdefwcs\fR, otherwise the first-defined user world system is used,
+else the physical system is used.
+
+If the MWCS is associated with an image then the "object" attribute of
+the physical system will return the name of the image or image section
+defined as the physical coordinate system for the MWCS.
+This is not necessarily the full image, e.g., in the case of
+a multidimensional image, the physical system might be any 2D plane of the
+image.  In the case of an event file image, the image name may include a
+filter or blocking factor.  References back to the raw data image based on
+MWCS physical coordinates will work so long as the raw image is opened
+using the name returned by the interface.  If the image is already open
+and was accessed by descriptor via MWCS, the descriptor may be retrieved
+by a \fImw_stati\fR call to fetch MW_IMDES.
+
+All MWCS coordinate systems have the standard attributes, with default values
+being supplied by the interface if not set by the application.  In particular,
+the logical and physical coordinate systems have attributes and may be
+treated as a special case of a world coordinate system by the application.
+
+.nh 3
+Axis Mapping
+
+    The coordinate transformation procedures (section 3.2) include support
+for a feature called \fIaxis mapping\fR, used to implement \fIdimensional
+reduction\fR.  A example of dimensional reduction occurs in IMIO, when
+an image section is used to specify a subraster of an image of dimension
+less than the full physical image.  For example, the section might specify
+a 1 dimensional line or column of a 2 or higher dimensional image, or a
+2 dimensional section of a 3 dimensional image.  When this occurs the
+applications sees a logical image of dimension equal to that of the image
+section, since logically an image section \fIis\fR an image.
+
+Dimensional reduction is implemented in MWCS by a transformation on the
+input and output coordinate vectors.  The internal MWCS coordinate system
+is unaffected by either dimensional reduction or axis mapping; axis mapping
+affects only the view of the WCS as seen by the application using the
+coordinate transformation procedures.
+
+For example, if the physical image is an image cube and we access the
+logical image section "[*,5,*]", an axis mapping may be set up which
+maps \fIphysical\fR axis 1 to logical axis 1, physical axis 2 to the
+constant 5, and physical axis 3 to logical axis 2.  The internal system
+remains 3 dimensional, but the application sees a 2 dimensional system.
+Upon input, the missing axis y=5 is added to the 2 dimensional input
+coordinate vectors, producing a 3 dimensional coordinate vector for
+internal use.  During output axis 2 is dropped and replaced by axis 3.
+
+The axis map is entered with \fImw_saxmap\fR and queried with \fImw_gaxmap\fR.
+Here, \fIaxno\fR is a vector, with \fIaxno[i]\fR specifying the logical axis
+to be mapped onto physical axis I.  If zero is specified the constant
+\fIaxval[i]\fR is used instead.  Axis mapping may be enabled or disabled
+with a call to \fImw_seti\fR.
+
+Axis mapping affects all of the coordinate transformation procedures,
+plus \fImw_translate\fR (since it defines a translation in terms of the
+logical system), and all of the coordinate system specification procedures
+having an "axis" parameter, e.g., \fImw_gwattrs\fR.
+Axis mapping is not used with those procedures which directly access or
+modify the physical or world systems (e.g., \fImw_slterm\fR or
+\fImw_swterm\fR) since full knowledge of the physical system is necessary
+for such operations.
+
+.nh 3
+Binding the Physical System
+
+    Recall that all coordinate systems are defined in terms of the physical
+system, and that the Lterm defines the mapping between the physical system
+and the logical system.  Transformations of the logical system leave the
+physical and world systems unaffected.  The only exception to this is the
+procedure \fImw_bindphys\fR, which binds the physical system to the current
+logical system, i.e., makes the current logical system the new physical system.
+This involves a transformation of the linear term (CD matrix) of each world
+system, since a world system is defined in terms of the physical system,
+and initialization of the Lterm to (normally) the identity matrix and zero
+translation vector.  This operation is irreversible, i.e., once
+\fImw_bindphys\fR is executed the original physical system is lost.
+
+In the case of an MWCS which is associated with an image opened with an image
+section, the new physical system is not strictly speaking the logical system,
+but the image matrix of the image being accessed, i.e,. the current image
+ignoring the image section.  Hence, following a call to \fImw_bindphys\fR,
+the Lterm will always describe the translation between the physical image
+matrix currently being accessed, and the logical system (image section).
+
+.nh 2
+Set/Stat Procedures
+
+    The MWCS status procedures, used to query or set the MWCS parameters,
+are as follows.
+
+.nf
+			mw_seti (mw, what, ival)
+		ival = mw_stati (mw, what)
+			mw_show (mw, outfd, what)
+.fi
+
+The currently defined interface parameters are the following.
+
+.nf
+	   Name       Type          Description
+
+	MW_AXMAP	b	enable or disable axis mapping
+	MW_IMDES	i	descriptor of associated image
+	MW_INTERP	i	interpolator type for sampled wcs
+	MW_NDIM		i	dimensionality of logical system
+	MW_NPHYSDIM	i	dimensionality of physical system
+	MW_NWCS		i	number of wcs defined
+	MW_WCS		i	currently active wcs
+.fi
+
+MW_NDIM may differ from MW_NPHYSDIM if dimensional reduction has been
+specified and axis mapping is enabled.  MW_NWCS returns the number of
+WCS currently defined; at least two WCS are always defined, i.e., the
+logical and physical systems (the world system will default to the
+physical system if not otherwise defined).  The index of the current
+default WCS is given by MW_WCS.  In the case of a sampled WCS, the
+interpolator type used by the coordinate transformation procedures is
+specified by MW_INTERP.
+
+.nh 2
+Utility Routines
+
+    The following routines are used internally within the interface to
+compile or evaluate transformations, and may be useful in applications
+code as well.
+
+.nf
+		 mw_invert[r|d] (o_ltm, n_ltm, ndim)
+		   mw_mmul[r|d] (ltm_1, ltm_2, ltm_out, ndim)
+		   mw_vmul[r|d] (ltm, ltv_in, ltv_out, ndim)
+		    mw_glt[r|d] (v1, v2, ltm, ltv, ndim)
+.fi
+
+These routines perform matrix inversion, multiplication of a matrix by another
+matrix, multiplication of a vector by a matrix, and general linear
+transformation (matrix multiply and addition of translation vector).
+
+.nh 2
+Datatypes and Precision
+
+    All floating point data is stored internally in MWCS using double
+precision.  Most of the interface procedures have both type real and type
+double versions, e.g., for entering Lterm or Wterm data.
+The single precision versions should be normally used unless double
+precision is required to represent the data.
+
+Although all floating point data is stored internally as type double,
+coordinate transformations performed at runtime may be carried out using
+either single or double precision computations, depending upon, e.g., whether
+\fImw_ctranr\fR or \fImw_ctrand\fR is called to perform the transformation.
+What happens is that when the transformation is compiled by \fImw_sctran\fR,
+two transformation descriptors are prepared, one for type real and the other
+for type double, with the appropriate descriptor being selected at runtime
+to carry out the transformation.  Hence the precision appropriate for the
+problem at hand can be employed without requiring that the worst case
+precision be used for all applications.
+
+.nh
+IMIO Interface to MWCS
+.nh 2
+Image Header Representation
+
+    When MWCS is used with image data, the encoded MWCS object is stored
+in the image header, and loaded into an MWCS descriptor when the image is
+accessed by an applications program.  The format in which the MWCS is
+stored in the image header depends upon the type of image.  If the image
+has a "flex-header" (as for QPOE and the new image structures) the MWCS
+is encoded in a machine independent binary format and stored in the image
+header as a variable length byte array.  This provides full generality
+and is the most efficient approach.
+
+For the older image formats which use a FITS header (OIF and STF) it is
+necessary to encode the MWCS as a series of FITS cards.  The proposed FITS
+WCS format, already in use for STF format images (with minor deviations
+from the standard), is used to represent the MWCS Wterm.  Additional FITS
+cards are necessary to represent the Lterm.  The (P,W) array for sampled
+WCS can also be represented in a FITS header, although this is awkward and
+inefficient if the number of samples is large.
+
+The FITS header keywords used to represent the Wterm, Lterm, and sampled
+WCS are the following.
+
+.nf
+	WCSDIM		WCS dimension (may differ from image)
+
+	CTYPEn		coordinate type
+	CRPIXn		reference pixel
+	CRVALn		world coords of reference pixel
+	CDi_j		CD matrix
+
+	CDELTn		CDi_i if CD matrix not used (input only)
+	CROTA2		rotation angle if CD matrix not used
+
+	LTVi		Lterm translation vector
+	LTMi_j		Lterm rotation matrix
+
+	WSVi_LEN	Number of sample points for axis I
+	WSVi_jjj	Sampled WCS array for axis I
+
+	WATi_jjj	WCS attributes for axis I
+.fi
+
+Contrary to MWCS convention, the WCS stored in a FITS format header
+defines the transformation from the logical system (image matrix)
+to the world system, rather than the physical system.  The MWCS Wterm
+is computed from the FITS representation by transforming the FITS WCS
+by the stored Lterm when the stored MWCS is loaded.
+
+The name format CDi_j varies slightly from the proposed FITS standard, but is
+backwards compatible with STF (and more readable than the FITS nomenclature).
+The keywords LTVECn, LTi_j, WSVi_LEN, WSVi_jjj, and WATi_jjj are peculiar
+to MWCS.  A sampled WCS is represented as a series of WSVi_jjj cards,
+wherein the sample points are stored as character strings, storing as
+many sample points as possible on each card, ignoring the card boundaries
+(i.e., a card may end in the middle of a number).  WCS attributes are likewise
+encoded as a series of WATi_jjj cards, giving the attributes for axis I as
+string data of the form "attribute = value", ignoring card boundaries.
+Multiple world coordinate systems (other than the physical and one world
+system) cannot be used with old format image headers.
+
+.nh 2
+Handling of the WCS by IMIO
+
+    When an image is opened by IMIO the image header is read, an MWCS
+descriptor is opened, and the stored MWCS is loaded into the MWCS descriptor
+from the image header.  If an image section has been opened the Lterm
+is then updated to reflect the additional linear transformation defined
+by the section.  The correct logical to physical or world transformation
+is then seen at the IMIO level, and will be propagated to a new image
+in a NEW_COPY image operation when the MWCS is copied to the new image.
+
+In the case of an image format which uses a FITS header, application of
+the section transform during an image open \fIdoes not\fR include updating
+of the FITS representation of the WCS.  There are two problems with doing
+so: all this editing of the FITS image of the header is inefficient
+unless absolutely necessary, and more seriously, if the image is opened
+READ_WRITE with an image section and the header is later updated, the
+stored WCS will be incorrect.  So, while the WCS as represented by the
+MWCS will always be correct, the FITS header parameters will reflect
+the WCS of the raw image ignoring the image section.  If it is necessary
+for some reason to update the FITS header in memory to reflect the image
+section, \fImw_saveim\fR may be called to perform the udpate.
+
+Propagation of the correct logical system in a NEW_COPY operation works
+because once the FITS header is copied, \fImw_saveim\fR is called to
+edit the header of the new image.
+
+.nh
+Implementation
+.nh 2
+Restrictions
+
+    Since there was not time to solve the general WCS problem with the MWCS
+interface, several restrictions were accepted for this version.  These are
+the following.
+.ls
+.ls o
+All WCS functions are built in (hard coded), hence the interface is not
+extensible at runtime and the only way to support new applications is
+through modification of the interface (by adding new function drivers).
+.le
+.ls o
+There is no support for modeling geometric distortions, except possibly
+in one dimension.
+.le
+.ls o
+There is no provision for storing more than one world coordinate system
+in FITS oriented image headers, although multiple WCS are supported internally
+by the interface, and are preserved and restored across \fImw_save\fR and
+\fImw_load\fR operations.
+.le
+.ls o
+Coordinate transforms involving dependent axes must includes all such axes
+explicitly in the transform.  Dependent axes are axes which are related,
+either by a rotation, or by a WCS function.  Operations which could subset
+dependent axis groups, and which are therefore disallowed, include setting
+up a transform with an AXES bitmap which excludes dependent axes, or more
+importantly, an image section involving dimensional reduction, where the
+axis to be removed is not independent.  This could happen, for example,
+if a two-dimensional image were rotated and one tried to open a
+one-dimensional section of the rotated image.
+.le
+.le
+
+All these problems can be solved given enough time, although the last problem
+mentioned becomes very complicated (perhaps intractable) when nonlinear world
+systems of dimension greater than three are involved.
+
+.nh 2
+Function Drivers
+
+    World coordinate systems are implemented in MWCS by providing something
+called a \fIfunction driver\fR for each function type, as specified by the
+\fIwtype\fR argument to \fImw_swtype\fR.  The \fIwtype\fR is the name of
+the function, and the name of the function driver.
+
+A function driver consists of the following procedures.  A given driver
+need not implement all driver procedures; procedures which are not used
+by a driver are set to NULL in the function driver table.
+
+.nf
+        operation                       syntax
+
+        FN_INIT                 wf_FCN_init (fc, dir)
+        FN_DESTROY           wf_FCN_destroy (fc)
+        FN_FWD                   wf_FCN_fwd (fc, pv, wv)
+        FN_INV                   wf_FCN_inv (fc, wv, pv)
+.fi
+
+where FCN is replaced by a 3 letter abbreviation for the function name,
+e.g., "smp" for the sampled WCS function, "tan" for the tangent plane
+projection etc.  This is only a suggested naming convention; the actual
+driver procedure names are arbitrary so long as name conflicts are
+avoided.
+
+The argument FC to each driver procedure is a pointer to the function call
+descriptor set up by \fImw_sctran\fR.  This consists of a number of standard
+fields followed by an area which is reserved for fields which are private
+to the function driver.  During compilation of a transformation, the function
+driver initialization procedure FN_INIT will be called to perform any function
+dependent initialization, e.g., processing of the attribute list for the
+axes assigned to the function, to input any function specific parameters.
+
+During runtime evaluation of a function call, FN_FWD will be called for a
+forward transformation (physical to world), and FN_INV for an inverse
+transformation (world to physical).  Note that the linear portion of the
+WCS, i.e., the CD matrix and all other linear terms except W (the CRVAL
+vector) are handled the same for all WCS functions, outside of the driver.
+Hence when the driver is called for a forward transformation, for example,
+the CD matrix and R vector (defining the reference point) will already
+have been applied to the input vector PV.
+
+To fully understand how function drivers are implemented it is probably
+simplest to study the existing drivers.
+
+.tp 40
+.sh
+Appendix A:  Interface Summary
+
+.nf
+		   mw = mw_open (bufptr|NULL, ndim)
+		 mw = mw_openim (im)
+		mw = mw_newcopy (mw)
+		       mw_close (mw)
+
+			mw_load (mw, bufptr)
+		  len = mw_save (mw, bufptr, buflen)
+	       mw_[load|save]im (mw, im)
+
+		 ct = mw_sctran (mw, system1, system2, axes)
+	  ndim = mw_gctran[r|d] (ct, ltm, ltv, axtype1, axtype2, maxdim)
+		      mw_ctfree (ct)
+
+	    x2 = mw_c1tran[r|d] (ct, x1)
+		 mw_v1tran[r|d] (ct, x1, x2, npts)
+		 mw_c2tran[r|d] (ct, x1,y1, x2,y2)
+		 mw_v2tran[r|d] (ct, x1,y1, x2,y2, npts)
+		  mw_ctran[r|d] (ct, p1, p2, ndim)
+		  mw_vtran[r|d] (ct, v1, v2, ndim, npts)
+
+	     mw_[s|g]lterm[r|d] (mw, ltm, ltv, ndim)
+	      mw_translate[r|d] (mw, ltv_1, ltm, ltv_2, ndim)
+		      mw_rotate (mw, theta, center, axes)
+		       mw_scale (mw, scale, axes)
+		       mw_shift (mw, shift, axes)
+
+		   mw_newsystem (mw, system, ndim)
+		 mw_[s|g]system (mw, system[, maxch])
+		  mw_[s|g]axmap (mw, axno, axval, ndim)
+		    mw_bindphys (mw)
+
+	     mw_[s|g]wterm[r|d] (mw, r, w, cd, ndim)
+		      mw_swtype (mw, axis, naxes, wtype, wattr)
+	     mw_[s|g]wsamp[r|d] (mw, axis, pv, wv, npts)
+		 mw_[s|g]wattrs (mw, axis, attribute, valstr[, maxch])
+
+		 mw_invert[r|d] (o_ltm, n_ltm, ndim)
+		   mw_mmul[r|d] (ltm_1, ltm_2, ltm_out, ndim)
+		   mw_vmul[r|d] (ltm, ltv_in, ltv_out, ndim)
+		    mw_glt[r|d] (v1, v2, ltm, ltv, ndim)
+
+			mw_seti (mw, what, ival)
+		ival = mw_stati (mw, what)
+			mw_show (mw, outfd, what)
+.fi
+.sp
+.endhelp
-- 
cgit 

