Subject: Comments on Draft WCS Standard

From: fred at manono.span (Fred Patt (301)286-4569)
Date: Thu, 18 Feb 93 12:59:38 EST

The following message was addressed to "fitsbits-request at nrao.edu",
rather than to "fitsbits at nrao.edu", and so it came to me, Don Wells,
the maintainer of the exploder, rather than being exploded as
intended.
 ------- Start of forwarded message ------- 
From: fred at manono.span (Fred Patt (301)286-4569) 
To: fitsbits-request at fits.CV.NRAO.EDU 
Subject: Comments on Draft WCS Standard 
Date: Thu, 18 Feb 93 11:12:32 EST

Dr. Greisen, 

	I was offered a chance by Richard White to review your draft World
Coordinate System standard for FITS, due to my interest in map projections
in general and the quad sphere in particular.  My recent involvement has
been with Earth science data sets but I was also involved with the COBE
mission for several years and am familiar with their skymap implementation.
I have advocated the adoption of the quad sphere by local Earth science 
projects and was responsible for the use of this scheme in the Advanced
Very High Resolution Radiometer (AVHRR) Pathfinder Land prototype
processing system here at Goddard.  I provided details of the projection
equations to the AVHRR group in the form of pseudocode, which has been
implemented in C by the developers.  The data system manager, Mary James of
Goddard Code 902, has given me encouragement to publish my work on this
project. 

	My point here is that an exact form of the projection was derived 
in the second of the two references (O'Neill and Laubscher, 1976) in terms 
of polar coordinates, and this is the form which was implemented by AVHRR.
The reference only derived the projection in the forward direction, from 
sky to cube coordinates, but I found that the derivation of the reverse 
projection was very straightforward.  I also found that the equations could 
be significantly simplified by direct use of cartesian coordinates, thereby
eliminating the polar angles as intermediaries and, as a side benefit, 
reducing the number of trig function evaluations required.  

	The information regarding the COBE implementation by Immanuel 
Freedman is completely accurate.  COBE has accepted the non-exact form of 
the transformation, but I suggest that the description of the exact form is 
at least as appropriate for your work.  The combination of the overall 
complexity of the projection equations and the errors in the projection 
described for COBE would seem to discourage most potential users.  

	I have attached a set of summary notes which I provided for the 
AVHRR group and the derivation of the reverse projection and cartesian form 
equations for the exact transformation (in computer-language syntax; it 
would be reasonable for me to put these into a more appropriate notation if 
desired).  I hope that this information may be of use to you.

						Fred Patt 

SUMMARY

	The quadrilaterlized spherical cube, or quad sphere for short, is 
an equal-area mapping and binning scheme for data collected on a spherical
surface (either Earth data or the celestial sphere).  It was first proposed
by Chan and O'Neill in 1975 for the Naval Environmental Prediction Research 
Facility (Reference 1).  
	
	There are two key elements to the quad sphere:

o  The mapping consists of projecting the sphere onto the faces of an 
inscribed cube using a curvilinear projection which preserves area.  The 
sphere is divided into six equal areas which correspond to the faces of the
cube.  The vertices of the cube correspond to the cartesian coordinates
defined by |x|=|y|=|z| on the unit sphere.  For an Earth projection, the 
cube is normally oriented with one face normal to the North Pole and one 
face centered on the Greenwich meridian (although any definition of pole
and meridian could be used).  The faces of the cube are divided into square
bins, where the number of bins along each edge is a power of 2, selected to
produce the desired bin size.  Thus the number of bins on each face is
2**(2*N), where N is the binning level, and the total number of bins is
6*2**(2*N).  For example, a level of 10 gives 1024x1024 bins on each face
and 6291456 (6*2**20) total bins, which are 23.605 square arcminutes 
(1.99737E-6 steradians) in size.

o  The bins are numbered serially, rather than being rastered as for an 
image.  The bin numbers are determined as follows.  The total number of
bits required for the bin numbers at level N is 2*N+3, where the 3 MSBs are
used for the face numbers and the remaining bits are used to number the
bins within each face.  The faces are numbered 0-5 with 0 being the North 
face, 1 through 4 being equatorial with 1 corresponding to Greenwich, and 5
being South.  Thus at level 10, face 0 has bin numbers 0-1048577, face 1 
has numbers 1048576-2097151, etc.  Within each face the bins are numbered 
serially from one corner (the convention is to start at the "lower left") 
to the opposite corner, with the ordering such that each pair of bits 
corresponds to a level of bin resolution.  This ordering in effect is a 
two-dimensional binary tree, which is referred to as the quad-tree  The 
conversion between bin numbers and coordinates is straightforward.  If 
4-byte integers are used for the bin numbers the maximum practical bin 
level is 14, which uses 31 of the 32 bits and results in a bin size of 
0.0922 square arcminutes (7.80223E-9 steradians).  The advantages of this 
numbering approach are stated below. 

	In principal the mapping and numbering schemes are separable; the 
projection onto the cube could be used with another bin numbering scheme, 
and the numbering scheme itself could be used with any arrangement of bins 
which can be partitioned as a set of square arrays.  Used together, they 
comprise a flexible and efficient system for storing map data.   

ADVANTAGES

	The quad sphere projection does not produce singularities at the 
poles or elsewhere, as do some other equal-area mapping schemes.  The 
distortion is moderate over the entire sphere, so that at no point are 
shapes distorted beyond recognition.  

	Individual faces can be used to display parts of the sphere without
remapping by converting from pixel order to raster order; coordinate grids
can be overlaid on the faces to assist in registration. 

	The bin numbering allows the resolution to be changed (by factors 
of two) simply by adding or deleting LSBs.  This is particularly useful if 
it is desired to increase the bin size, for example to compare maps with 
different resolutions.  This is performed simply by dividing the bin
numbers by four for each factor of two increase in bin size.  

	The bin numbers tend be closer together in sequence for neighboring 
bins.  

	The quadtree bin numbering allows maps to be stored as FITS
BINTABLES with no wasted storage needed for blank areas (unfilled bins).
Only bins containing valid data need be included in a file, as the bin
numbers themselves identify the locations of the bins. 
	
	The serial numbering frees up one array dimension for another use
(e.g., time, a third coordinate, multichannel data, etc.). 

DERIVATIONS OF REVERSE AND CARTESIAN TRANSFORMATIONS 

The equations for the polar form of the transformation from the sphere to 
the cube are given in reference 2 as equations (3-21) and (3-38).  I have 
taken this approach further in two areas:  1) inversion of these equations 
to derive the cube-to-sphere transformation; and 2) simplification of these 
equations to perform the transformation directly in terms of cartesian 
coordinates.  Both of these are straightforward as shown below.

1.  DERIVATION OF THE CUBE-TO-SPHERE TRANSFORMATION

The exact form of the sphere-to-cube tranformation in reference 2 is 
derived in terms of the co-elevation and azimuth angles (theta, phi) on the 
sphere and an equivalent set of angles (mu,nu) on the cube.  The equal-area
tranformation is performed using the following equations:

TAN(mu) = (12/pi)*(theta+ACOS(SIN(theta)*SIN(pi/4))-pi/2) 	(3-21)

TAN(nu) = SQRT((SEC(mu))**2*(1-COS(theta))/(1-COS(ATAN(SEC(theta)))))  (3-38)

The inversion of these equations is as follows.  Using the trig identity 
ASIN(x)+ACOS(x)=pi/2, and the knowledge that SIN(pi/4) = 1/SQRT(2), 
equation (3-21) can be written as

TAN(mu) = (12/pi)*(theta-ASIN(SIN(theta)/SQRT(2)))		(3-21a)

and then rearranged as

theta - (pi/12)*TAN(mu) = ASIN(SIN(theta)/SQRT(2))

Taking the sine of both sides gives

SIN(theta)*COS((pi/12)*TAN(mu)) - COS(theta)*SIN((pi/12)*TAN(mu))
	 = SIN(theta)/SQRT(2)

Dividing through by COS(theta) and solving for TAN(theta) yields the 
desired inverse relationship:

TAN(theta) = SIN((pi/12)*TAN(mu))/(COS((pi/12)*TAN(mu))-1/SQRT(2)) 	(1)

The inversion of equation (3-38) is performed by squaring both sides and 
solving for COS(phi):

COS(phi) = 1 - (TAN(nu))**2*(1-COS(ATAN(SEC(theta))))/(SEC(mu)**2	(2)

As in the original derivation, this equation used the prior determination 
of theta as a function of mu.  

2.  EXPRESSION OF THE TRANSFORMATIONS IN CARTESIAN COORDINATES

Both sets of transformation equations can be expressed in cartesian 
coordinates (three-dimensional coordinates of the unit vector on the sphere 
and two-dimensional coordinates on the cube face).  Assume the unit vector 
has coordinates (q,r,s) and the coordinates on the cube face are (x,y),
where q is normal to the plane of the cube face and (r,s) cooresponds to
(x,y).  Then the functions of angles can be expressed as follows: 

	TAN(theta) = s/r
	SIN(theta) = s/SQRT(r*r+s*s)
	COS(phi)   = q
	TAN(mu)    = y/x
	TAN(nu)    = SQRT(x*x+y*y)

Substituting these into the transformation equations results in the 
following revised forms for the equations.

Equation (3-21a):

	y/x = (12/pi)*(ATAN(s/r)- ASIN(s/SQRT(2*r*r+2*s*s)))	(3)

Equation (3-38) can first be simplified by the following substitution:

	COS(ATAN(SEC(theta))) = 1/SQRT(1+(SEC(theta))**2)

Substituting this and the appropriate cartesian coordinate expressions 
gives: 

	x*x+y*y = (1+(y/x)**2)*(1-q)/(1-1/SQRT(2+(s/r)**2))

which can be re-written to eliminate y:

	x = sqrt((1-q)/(1-1/SQRT(2+(s/r)**2)))			(4)

Note that there is no ambiquity in the sign of the square root since the 
derivation of (3-38) assumes that x is positive.  This and equation (3) for 
y/x above give the solution for both y and x, noting also that (3-21) 
assumes that x is not less than y.  

Equation (1) above can be rewritten also by direct substitution:

	s/r = SIN((pi/12)*(y/x))/(COS((pi/12)*(y/x))-1/SQRT(2))	    (5)

Equation (2) above can be rewritten by solving (4) for q:

	q = 1 - x*x*(1-1/SQRT(2+(s/r)**2)))			(6)

Thus equations (3) through (6) can be used to convert directly from 
coordinates on the sphere to the cube and back without the need to convert 
to angles as intermediate variables, and the number of trigonometric 
function evaluations required is greatly reduced. 

REFERENCES

1.  Chan, F.K. and O'Neill (1975), Feasibility Study of a Quadrilateralized
Spherical Cube Earth Data Base, Computer Sciences Corporation, EPRF
Technical Report 2-75 (CSC).  Prepared for the Environmental Prediction
Research Facility, Monterey, California. 

2.  O'Neill, I.M. and Laubscher, R.E. (1976), Extended Studies of a
Quadrilateralized Spherical Cube Earth Data Base, Computer Sciences
Corporation, NEPRF Technical Report 3-76 (CSC). Prepared for the Naval
Environmental Prediction Research Facility, Monterey, California. 
------- End of forwarded message -------