Coordinate.h

Go to the documentation of this file.

//# Coordinate.h: Interface for converting between world and pixel coordinates
//# Copyright (C) 1997,1999,2000,2001,2002,2003,2004
//# Associated Universities, Inc. Washington DC, USA.
//#
//# 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 Massachusetts Ave, Cambridge, MA 02139, USA.
//#
//# Correspondence concerning AIPS++ should be addressed as follows:
//#        Internet email: aips2-request@nrao.edu.
//#        Postal address: AIPS++ Project Office
//#                        National Radio Astronomy Observatory
//#                        520 Edgemont Road
//#                        Charlottesville, VA 22903-2475 USA
//#
//#
//# $Id$


#ifndef COORDINATES_COORDINATE_H
#define COORDINATES_COORDINATE_H

#include <casacore/casa/aips.h>
#include <casacore/casa/BasicSL/String.h>
#include <casacore/casa/Arrays/Vector.h>
#include <wcslib/wcs.h>

namespace casacore { //# NAMESPACE CASACORE - BEGIN


template<class T> class Quantum;
template<class T> class Matrix;
class IPosition;
class RecordInterface;
class Projection;

// <summary>
// Interface for converting between world and pixel coordinates.
// </summary>

// <use visibility=export>

// <reviewed reviewer="Peter Barnes" date="1999/12/24">
// </reviewed>

// <prerequisite>
//   <li> Knowledge of astronomical coordinate conversions in general. Probably the
//        best documents are the papers by Mark Calabretta and Eric Greisen.
//        The initial draft from 1996 can be found at
//        http://www.atnf.csiro.au/~mcalabre.  It is this draft that the
//        Coordinate classes are based upon.  Since then, this paper has evolved
//        into three which can be found at the above address, and will be published in the
//        Astronomy and Astrophysics Supplement Series (probably in 2000).
//        The design has changed since the initial draft.  When these papers
//        are finalized, and the IAU has ratified the new standards, WCSLIB
//        (Mark Calabretta's implementation of these conventions) will be
//        revised for the new designs.  At that time, the Coordinate classes
//        may also be revised.
//   <li> Generic Casacore classes; especially those in the 
//        <linkto module=Arrays>Arrays</linkto> module.
//   <li> Perhaps some of the information in the
//        <linkto module=Measures>Measures</linkto> module.
// </prerequisite>
//
// <synopsis>
// The Coordinate class defines the generic interface whereby a pixel position
// is converted to a world (sky, frequency, stokes, ...) position and vice
// versa. The pixel and world coordinates are in general
// multi-dimensional values. In general there need not be the same number of
// pixel and world axes, although this will normally be the case.
//
// The fundamental model is that a pixel is first converted into a relative
// physical coordinate by:
// <ol>
//    <li> Subtracting a reference pixel value from the pixel location; then
//    <li> Multiplying this offset by a general transformation matrix (usually
//         to account for rotation, but any matrix is allowed); then
//    <li> Multiplying this product by an increment in physical units.
// </ol>
// After this linear stage, the final coordinate value is computed from this
// relative physical unit and a reference value, and possibly some other 
// parameters. In the case of a sky position, these latter include at least the
// projection type. In the case of a purely linear coordinate, the reference value
// is merely added to the relative physical coordinate. The interface also
// allows the axes to be assigned names (reasonable defaults will be selected),
// and for physical units.
// 
// Both absolute and relative coordinates are supported.  The main
// interface supports conversion between absolute pixel
// and absolute world coordinate.  There are then functions to
// convert absolute coordinates to relative and vice versa.
// A relative pixel coordinate is defined according to 
// 
// relative = absolute - reference
//
// A relative world coordinate is similar, although there may
// be deviations from this formula (e.g. for DirectionCoordinate
// a cos(latitude) term is incorporated and for StokesCoordinate
// relative world coordinates  are defined to be the same as
// absolute world coordinates.
//
// </synopsis>
//
// <note role=caution>
// All absolute pixels coordinates are zero relative.
// </note>
//
// <example>
// This is a base class so there is no direct example, but
// see the example in <linkto module=Coordinates>Coordinates.h</linkto>
// for use of the derived classes.
// </example>
//
// <motivation>
// Encapsulate the common interface to coordinate conversion so that it may
// be used polymorphically.
// </motivation>
//
// <thrown>
//   <li>  AipsError
//   <li>  AllocError
// </thrown>
//
// <todo asof="1997/1/13">
//   <li> Perhaps common FITS related interfaces should go in this class.
// </todo>
//

class Coordinate
{
public:
    // This enum lists the types of the derived classes. It is primarly used
    // in the CoordinateSystem class.
    enum Type { 
        // Linear axes.
        LINEAR, 
        // A direction. Usually RA/DEC.
        DIRECTION, 
        // A spectral axis.
        SPECTRAL, 
        // A Stokes axis.
        STOKES, 
        // A one-dimensional Cooordinate system, usually created from a table
        // although it can also be purely linear.
        TABULAR,
        // to mark DATA and ERROR values
        QUALITY,
        // A CoordinateSystem (a collection of Coordinates).
        COORDSYS };

    // This enum is used for formatting world values into Strings
    enum formatType {
       // Default; formatter decides
       DEFAULT,
       // Scientific format (e.g. -1.2397E+03)
       SCIENTIFIC,
       // Fixed floating format (e.g. 12.134)
       FIXED,
       // Either scientific or floating point, auto-selected by the C++
       // STL formatting routines.  May not be available for all Coordinate
       // types.
       MIXED,
       // HHH:MM:SS.SSS style formatting
       TIME };

    // Destructor.  Needs to be public so the user can delete Coordinate* objects
    virtual ~Coordinate();

    // List the type of this Coordinate object. 
    // <group>
    virtual Type type() const = 0;
    virtual String showType() const = 0;
    static String typeToString (Coordinate::Type type);
    // </group>

    // How many world/pixel axes are there in this Coordinate? While the number
    // of world and pixel axes will generally be the same, it is not a 
    // requirement. For example, in CoordinateSystem you could remove a pixel
    // axis and leave the corresponding world axis. Also, if we ever implement
    // a "SlicedCoordinate" class then there would be more world than pixel
    // coordinates (the pixel coordinate would be a pixel number along the slice,
    // whereas the world axes would continue to be RA/DEC).
    // <group>
    virtual uInt nPixelAxes() const = 0;
    virtual uInt nWorldAxes() const = 0;
    // </group>

    // Convert an absolute pixel position to an absolute world position or vice 
    // versa. Returns True
    // if the conversion succeeds, otherwise it returns False and method
    // errorMessage contains an error message. The input vector must be of length
    // <src>nPixelAxes</src> or <src>nWorldAxes</src>.  The output vector
    // is resized appropriately.
    // if <src>useConversionFrame</src>, if the coordinate has a conversion layer frame
    // (such as can be present in spectral and direction coordinates), it
    // is used. Else, the native frame is used for the conversion.
    // <group>
    virtual Bool toWorld(Vector<Double> &world, 
                         const Vector<Double> &pixel, Bool useConversionFrame=True) const = 0;
    virtual Bool toPixel(Vector<Double> &pixel, 
                         const Vector<Double> &world) const = 0;
    // </group>

    // Mixed absolute pixel/world coordinate conversion.
    // worldIn and worldAxes are vectors of length <src>nWorldAxes</src>.
    // <src>pixelIn</src> and <src>pixelAxes</src> are of length <src>nPixelAxes</src>.
    // <src>worldAxes(i) = True</src> specifies you have given a world
    // value in <src>worldIn(i)</src> to convert to pixel.
    // <src>pixelAxes(i)=True</src> specifies you have given a pixel 
    // value in <src>pixelIn(i)</src> to convert to world.
    // You cannot specify the same axis via <src>worldAxes</src>
    // and <src>pixelAxes</src>.
    // Values in <src>pixelIn</src> are converted to world and
    // put into <src>worldOut</src> in the appropriate world axis
    // location.  Values in <src>worldIn</src> are copied to
    // <src>worldOut</src>.
    // Values in <src>worldIn</src> are converted to pixel and
    // put into <src>pixelOut</src> in the appropriate pixel axis
    // location.  Values in <src>pixelIn</src> are copied to
    // <src>pixelOut</src>.
    // <src>worldMin</src> and <src>worldMax</src> specify the range of the world
    // coordinate (in the world axis units of that world axis
    // in the CoordinateSystem) being solved for in a mixed calculation
    // for each world axis. They are only actually needed for DirectionCoordinates
    // and for all other Coordinates the relevant elements   
    // can be undefined.   If you don't know, use -180 to 180
    // degrees for longitude, and -90 to 90 for latitude.
    // Removed axes are handled (for example, a removed pixel
    // axis with remaining corresponding world axis will
    // correctly be converted to world using the replacement
    // value).
    // Returns True if the conversion succeeds, otherwise it returns False and
    // <src>errorMessage()</src> contains an error message. The output vectors
    // are resized.
    virtual Bool toMix(Vector<Double>& worldOut,
                       Vector<Double>& pixelOut,
                       const Vector<Double>& worldIn,
                       const Vector<Double>& pixelIn,
                       const Vector<Bool>& worldAxes,   
                       const Vector<Bool>& pixelAxes,
                       const Vector<Double>& worldMin,
                       const Vector<Double>& worldMax) const;

    // Set the world min and max ranges, for use in function <src>toMix</src>, for 
    // a lattice of the given shape for this coordinate. The default implementation
    // here sets the range for pixels dangling 25% off the image.
    // Returns False if fails with a reason  in <src>errorMessage()</src>.
    // setDefaultWorldMixRanges sets the range for each axis to +/-1e99
    // The ranges remain zero length vectors until you explicitly
    // initialize them.
    // <group>
    virtual Bool setWorldMixRanges (const IPosition& shape);
    virtual void setDefaultWorldMixRanges ();
    Vector<Double> worldMixMin () const {return worldMin_p;};
    Vector<Double> worldMixMax () const {return worldMax_p;};
    //</group>


    // Batch up a lot of transformations. The first (most rapidly varying) axis
    // of the matrices contain the coordinates. Returns False if any conversion
    // failed  and  <src>errorMessage()</src> will hold a message.
    // The <src>failures</src> array (True for fail, False for success)
    // is the length of the number of conversions and
    // holds an error status for each conversion.  The default
    // implementation is provided that works with the "single" version of
    // <src>toWorld</src> and <src>toPixel</src>, but for maximum efficiency these should be
    // overridden.
    // <group>
    virtual Bool toWorldMany(Matrix<Double>& world, 
                             const Matrix<Double>& pixel, 
                             Vector<Bool>& failures) const;
    virtual Bool toPixelMany(Matrix<Double>& pixel, 
                             const Matrix<Double>& world, 
                             Vector<Bool>& failures) const;
    // </group>

    // Make absolute coordinates relative and vice-versa (with
    // respect to the reference value).
    // Vectors must be length <src>nPixelAxes()</src> or
    // <src>nWorldAxes()</src> or memory access errors will occur
    // <group>
    virtual void makePixelRelative (Vector<Double>& pixel) const;
    virtual void makePixelAbsolute (Vector<Double>& pixel) const;
    virtual void makeWorldRelative (Vector<Double>& world) const;
    virtual void makeWorldAbsolute (Vector<Double>& world) const;
    // </group>

    // Make absolute coordinates relative and vice versa with respect
    // to the given reference value.  Add the other functions in this grouping
    // as needed. Vectors must be length <src>nPixelAxes()</src> or
    // <src>nWorldAxes()</src> or memory access errors will occur
    // <group>
    virtual void makeWorldAbsoluteRef (Vector<Double>& world,
                                       const Vector<Double>& refVal) const;
    // </group>


    // Batch up a lot of absolute/relative transformations. 
    // Parameters as above  for 
    // <src>toWorldMany</src> and <src>toPixelMany</src>
    // <group>
    virtual void makePixelRelativeMany (Matrix<Double>& pixel) const;
    virtual void makePixelAbsoluteMany (Matrix<Double>& pixel) const;
    virtual void makeWorldRelativeMany (Matrix<Double>& world) const;
    virtual void makeWorldAbsoluteMany (Matrix<Double>& world) const;
    // </group>


    // Return the requested attributed.
    // <group>
    virtual Vector<String> worldAxisNames() const = 0;
    virtual Vector<Double> referencePixel() const = 0;
    virtual Matrix<Double> linearTransform() const = 0;
    virtual Vector<Double> increment() const = 0;
    virtual Vector<Double> referenceValue() const = 0;
    virtual Vector<String> worldAxisUnits() const = 0;
    // </group>

    // Set the requested attribute.  Note that these just
    // change the internal values, they do not cause any recomputation.
    // <group>
    virtual Bool setWorldAxisNames(const Vector<String> &names) = 0;
    virtual Bool setReferencePixel(const Vector<Double> &refPix) = 0;
    virtual Bool setLinearTransform(const Matrix<Double> &xform) = 0;
    virtual Bool setIncrement(const Vector<Double> &inc)  = 0;
    virtual Bool setReferenceValue(const Vector<Double> &refval)  = 0;
    // </group>

    // Change the units. Adjust the increment and
    // reference value by the ratio of the old and new units. This implies that
    // the units must be known <linkto class=Unit>Unit</linkto> strings, and that
    // they must be compatible, e.g. they can't change from time to length.
    //
    // A default implementation is available which does everything except set
    // the units vector, which must be done in the derived class.
    virtual Bool setWorldAxisUnits(const Vector<String> &units) = 0;

    // Find the Coordinate for when we Fourier Transform ourselves.  This pointer 
    // must be deleted by the caller. Axes specifies which axes of the Coordinate
    // you wish to transform.   Shape specifies the shape of the image
    // associated with all the axes of the Coordinate. Currently the
    // output reference pixel is always shape/2.
    virtual Coordinate* makeFourierCoordinate (const Vector<Bool>& axes,
                                               const Vector<Int>& shape) const;

    // If the last conversion to world or pixel coordinates resulted in an
    // error, report that error. If the last conversion succeeded, it is
    // undefined what this will return (it might well contain the last error
    // message).
    const String& errorMessage() const;

    // Comparison to fractional tolerance (for floating point values). 
    // Don't compare on specified axes in Coordinate. If the comparison
    // returns False, <src>errorMessage()</src> contains a message.
    // <group>
    virtual Bool near(const Coordinate& other, 
                      Double tol=1.0e-6) const = 0;
    virtual Bool near(const Coordinate& other, 
                      const Vector<Int>& excludeAxes,
                      Double tol=1.0e-6) const = 0;
    // </group>


    // Provide a common interface to getting formatted representations of
    // coordinate values.    Different derived Coordinate types are formatted
    // in different ways.  For example, an RA/DEC  DirectionCoordinate
    // uses an HMS.SS/DMS.SS representation. A Galactic Lat/Long DirectionCoordinate
    // uses floating format in degrees.  Other derived Coordinates are formatted with 
    // scientific format or floating format. The derived class format functions
    // provide this functionality.   
    // 
    // You may specify the format with the format argument and a value
    // from the enum <src>Coordinate::formatType</src>. If you give it the value 
    // <src>Coordinate::DEFAULT</src> then a sensible default is used.
    //
    // A mechanism for specifying the precision number of significant digits after 
    // decimal point is provided.  You can specify the precision directly when 
    // calling format if it is unambiguous how the derived Coordinate is 
    // going to be formatted.  For example, a LinearCoordinate is always formatted with 
    // scientific format.  However, if you are using these classes polymorphically, you 
    // don't want to have to know this and some derived Coordinates may be formatted
    // in multiple ways (such as the DirectionCoordinate examples above).
    // Therefore, the function getPrecision enables 
    // you to set default precisions for the different styles of formatting 
    // used variously in the base and derived classes.   This function chooses the 
    // precision from these default values, according to the type of derived 
    // Coordinate that your object is and what value for format that
    // you give (refer to the derived classes for details on this).
    // 
    // Some derived classes will format differently depending upon whether
    // you want to format an absolute or offset world value input via 
    // absolute (e.g. DirectionCoordinates).
    // 
    // The provided <src>worldValue</src> must be in the native units
    // of the Coordinate.  It may be an absolute (<src>isAbsolute=True</src>)
    // or relative (<src>isAbsolute=False</src>) value.   You may choose to
    // format the world value as absolute (<src>showAsAbsolute=True</src>) or
    // relative (<src>showAsAbsolute=False</src>).  <src>axis</src>
    // specifies which axis of the Coordinate this value belongs to.
    //
    // <src>units</src> specifies the units in which the input world value
    // will be formatted.  
    // If <src>units</src> is empty, the native unit for the given axis 
    // is used.
    //
    // Some derived classes will format in units different from the 
    // native unit of the Coordinate. The units of
    // the formatted number are returned in <src>units</src>.
    // If the <src>units</src> string is provided, the unit must be
    // consistent with the native unit of the coordinate.  The input
    // world value will be converted to this unit.
    //
    // You can also use the Quantum interface.  The units of the Quantum 
    // can then be anything  consistent with the Coordinate.
    // 
    // The default implementation here is to format only
    // with scientific or fixed formats. If precision is negative, a 
    // the default precision is used.
    //
    //<group>
    virtual void getPrecision(Int &precision,
                              Coordinate::formatType& format,
                              Bool showAsAbsolute,
                              Int defPrecScientific,
                              Int defPrecFixed,
                              Int defPrecTime) const;
    virtual String format(
        String& units,
        Coordinate::formatType format,
        Double worldValue,
        uInt axis,
        Bool isAbsolute=True,
        Bool showAsAbsolute=True,
        Int precision=-1,
        Bool usePrecForMixed=False
    ) const;

    String formatQuantity(String& units,
                          Coordinate::formatType format, 
                          const Quantum<Double>& worldValue, 
                          uInt axis, 
                          Bool isAbsolute=True,
                          Bool showAsAbsolute=True,
                          Int precision=-1);
    //</group>

    // Used for persistence. Derived classes will have similar static
    // restore methods. It will typically only return False if fieldName
    // has already been defined.
    virtual Bool save(RecordInterface &container,
                    const String &fieldName) const = 0;

    // Make a copy of ourself. This pointer has been allocated with
    // <src>new</src> and must be deleted by the caller.
    virtual Coordinate *clone() const = 0;

    // Comparison only made for specified axes in this and other Coordinate 
    // The default implementation should be ok for all Coordinate types
    // except Stokes and Quality...
    virtual Bool doNearPixel (const Coordinate& other, 
                              const Vector<Bool>&  thisAxes,
                              const Vector<Bool>& otherAxes,
                              Double tol=1.0e-6) const;

    // return the result of rotating the coordinate clockwise through the specified angle.
    // Rotation occurs about the reference pixel.
    // Coordinate must have exactly two pixel axes. The return type is the same
    // as the input type. It is the caller's responsibility to delete the returned pointer
    // when done with it to prevent a memory leak.
    // This method ultimately just changes the input coordinate's linear transform matrix.
    virtual Coordinate* rotate(const Quantum<Double>& angle) const;

    // Call wcsset on the wcs structure
    static void set_wcs (::wcsprm& wcs);

    // Call wcsini on the wcs structure
    static void init_wcs (::wcsprm& wcs, int naxis);

    // Call wcssub on the src/dst pair
    static void sub_wcs(const ::wcsprm &src, int &nsub, int axes[], ::wcsprm &dst);

    // Call wcssub on the src/dst pair with null nsub/axes
    static void copy_wcs(const ::wcsprm &src, ::wcsprm &dst);

protected:
    // Default constructor. Make an empty coordinate.  Used by derived classes.
    Coordinate();

    // Copy constructor (copy semantics)
    Coordinate(const Coordinate& other);

    // Assignment (copy semantics) 
    Coordinate& operator=(const Coordinate& other);

    // Set error message
    void set_error(const String &errorMsg) const;

    //
    Bool find_scale_factor(String &error, Vector<Double> &factor, 
                           const Vector<String> &units, 
                           const Vector<String> &oldUnits);


    // Tries to find a canonical unit for input unit (e.g.  GHz -> Hz), and
    // tells you the output name and unit for the Fourier coordinate 
    // pairing with the canonical unit
    void fourierUnits (String& nameOut, String& unitOut, String& unitInCanon,
                       Coordinate::Type type, Int axis, 
                       const String& unitIn, 
                       const String& nameIn) const;

   // Functions to interconvert pixel<->world via wcs.  These functions are called 
   // explicitly by the to{world,Pixel} functions in the appropriate wcs-based derived
   // classes. 
   // <group>
   Bool toWorldWCS (Vector<Double> &world, const Vector<Double> &pixel, wcsprm& wcs) const;
   Bool toPixelWCS(Vector<Double> &pixel,  const Vector<Double> &world, wcsprm& wcs) const;
   Bool toWorldManyWCS (Matrix<Double>& world, const Matrix<Double>& pixel,
                        Vector<Bool>& failures, wcsprm& wcs) const;
   Bool toPixelManyWCS (Matrix<Double>& pixel, const Matrix<Double>& world,
                        Vector<Bool>& failures, wcsprm& wcs) const;

   // Functions for handling conversion between the current units and
   // the wcs units. These are called explicitly by the appropriate 
   // derived class.
   // <src>convertFrom</src>
   // <group>
   void toCurrentMany (Matrix<Double>& world, const Vector<Double>& toCurrentFactors) const;
   void fromCurrentMany(Matrix<Double>& world, const Vector<Double>& toCurrentFactors) const;
   // </group>


   // Functions for handling conversion between the current reference frame 
   // and the native one. The default implementations do nothing.  They
   // should be over-ridden in the derived classes.
   // <group>
   virtual void convertTo (Vector<Double>&) const
     {}
   virtual void convertFrom (Vector<Double>&) const
     {}
   // </group>

   // Functions for handling conversion between the current reference frame 
   // and the native one for many conversions.  These functions just
   // call the virtual functions for single conversions.
   // <group>
   void convertToMany (Matrix<Double>& world) const;
   void convertFromMany (Matrix<Double>& world) const;
   // </group>

   // Interconvert between wcs PC cards and Matrix xForm format 
   void pcToXform (Matrix<Double>& xForm, const wcsprm& wcs) const;
   void xFormToPC (wcsprm& wcs, const Matrix<Double>& xForm) const;
   // </group>

    // toMix ranges.  Should be set by derived class.
    Vector<Double> worldMin_p, worldMax_p;

private:
    mutable String error_p;

    // Check format type
    void checkFormat(Coordinate::formatType& format,         
                     const Bool absolute) const;

    void makeWorldAbsRelMany (Matrix<Double>& value, Bool toAbs) const; 
    void makePixelAbsRelMany (Matrix<Double>& value, Bool toAbs) const; 


};

//###### Inlines

inline const String& Coordinate::errorMessage() const
{
    return error_p;
}

} //# NAMESPACE CASACORE - END

#endif