CoordinateSystem.h

Go to the documentation of this file.
00001 //# CoordinateSystem.h: Interconvert pixel and image coordinates.
00002 //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003,2004
00003 //# Associated Universities, Inc. Washington DC, USA.
00004 //#
00005 //# This library is free software; you can redistribute it and/or modify it
00006 //# under the terms of the GNU Library General Public License as published by
00007 //# the Free Software Foundation; either version 2 of the License, or (at your
00008 //# option) any later version.
00009 //#
00010 //# This library is distributed in the hope that it will be useful, but WITHOUT
00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00012 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00013 //# License for more details.
00014 //#
00015 //# You should have received a copy of the GNU Library General Public License
00016 //# along with this library; if not, write to the Free Software Foundation,
00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00018 //#
00019 //# Correspondence concerning AIPS++ should be addressed as follows:
00020 //#        Internet email: aips2-request@nrao.edu.
00021 //#        Postal address: AIPS++ Project Office
00022 //#                        National Radio Astronomy Observatory
00023 //#                        520 Edgemont Road
00024 //#                        Charlottesville, VA 22903-2475 USA
00025 //#
00026 //#
00027 //# $Id$
00028 
00029 
00030 #ifndef COORDINATES_COORDINATESYSTEM_H
00031 #define COORDINATES_COORDINATESYSTEM_H
00032 
00033 #include <casacore/casa/aips.h>
00034 #include <casacore/coordinates/Coordinates/Coordinate.h>
00035 #include <casacore/measures/Measures/MDirection.h>
00036 #include <casacore/measures/Measures/MFrequency.h>
00037 #include <casacore/coordinates/Coordinates/ObsInfo.h>
00038 #include <casacore/casa/Containers/Block.h>
00039 #include <casacore/measures/Measures/MDoppler.h>
00040 
00041 namespace casacore { //# NAMESPACE CASACORE - BEGIN
00042 
00043 template<class T> class Matrix;
00044 class DirectionCoordinate;
00045 class LinearCoordinate;
00046 class SpectralCoordinate;
00047 class StokesCoordinate;
00048 class QualityCoordinate;
00049 class TabularCoordinate;
00050 class IPosition;
00051 class LogIO;
00052 
00053 
00054 // <summary>
00055 // Interconvert pixel and world coordinates.
00056 // </summary>
00057 
00058 // <use visibility=export>
00059 
00060 // <reviewed reviewer="Peter Barnes" date="1999/12/24" tests="tCoordinateSystem">
00061 // </reviewed>
00062 //
00063 // <prerequisite>
00064 //   <li> <linkto class=Coordinate>Coordinate</linkto>
00065 // </prerequisite>
00066 
00067 // <synopsis>
00068 // CoordinateSystem is the normal interface to coordinate systems,
00069 // typically attached to an 
00070 // <linkto class=ImageInterface>ImageInterface</linkto>, however the
00071 // coordinate system can be manipulated on its own. CoordinateSystem
00072 // is in turn composed from various classes derived from the base class
00073 // <linkto class=Coordinate>Coordinate</linkto>.
00074 // <p>
00075 // The fundamental operations available to the user of a 
00076 // CoordinateSystem are:
00077 // <ol>
00078 //   <li> Transform a world (physical) coordinate to a pixel coordinate 
00079 //        or vice versa via the methods toWorld and toPixel.
00080 //   <li> Compose a CoordinateSystem from one or more independent groups,
00081 //        typically the sky-plane transformation will be one group, and the
00082 //        spectral axis will be another group. Each group consists of a linear
00083 //        transformation (in FITS terms, apply <src>CRPIX, PC, CDELT</src>)
00084 //        to turn the pixel coordinates into relative world coordinates, 
00085 //        followed by a (possibly) nonlinear projection to world coordinates 
00086 //        (i.e. apply <src>CTYPE and CRVAL</src>), typically a sky projection
00087 //        or a frequency to velocity conversion. Note that an arbitrary rotation
00088 //        or linear transformation can be applied by changing the
00089 //        matrix.
00090 //   <li> Transpose the world and/or pixel axes.
00091 //   <li> One or more pixel or world axes may be removed. You are encouraged to
00092 //        leave all the world axes if you remove a pixel axis.
00093 //        Removing a world axis also removes the corresponding pixel axis.
00094 //   <li> Calculate the CoordinateSystem that results from a subimage
00095 //        operation.
00096 // </ol>
00097 //
00098 // Note that all the knowledge to do with removing and transposing axes is
00099 // maintained by the CoordinateSystem.  The individual Coordinates, of which it
00100 // is made, know nothing about this.
00101 // <p>
00102 // Although the CoordinateSystem exists in the absence of an image, the usual
00103 // place you will find one is attached to an object derived from ImageInterface 
00104 // such as PagedImage. When you do so, the physical (or pixel) axes in the image
00105 // map one to one with the pixel axes contained in the CoordinateSystem.
00106 // It cannot be any other way as when you create a PagedImage, it is checked
00107 // that there are equal numbers of image and CoordinateSystem pixel axes.
00108 // It is up to the creator of the PagedImage to make sure that they are
00109 // in the correct order.
00110 // <p>
00111 // However, the CoordinateSystem may have more world axes than pixel axes
00112 // because it is possible to remove a pixel axis but not its associated
00113 // world axis (for example for a moment image).   Now, if you use
00114 // the CoordinateSystem functions 
00115 // referencePixel and referenceValue, you will find the vector of reference
00116 // values will have more values than the vector of reference pixels,
00117 // if a pixel axis has been removed but not the world axis.  You 
00118 // must use the ancilliary functions provided
00119 // to find out what is where.   
00120 // <p>
00121 // Let's consider an example where a CoordinateSystem consisted of
00122 // a DirectionCoordinate and a SpectralCoordinate.  Let us say that
00123 // the first two pixel axes of the image associate (roughly of course
00124 // because lines of constant RA and DEC are not parallel with
00125 // the pixel coordinates) with the DirectionCoordinate (RA and DEC say) 
00126 // and the third pixel axis is the SpectralCoordinate.
00127 // Now imagine we collapse the image along the second pixel axis (roughly,
00128 // the DEC axis).  For the output image, we remove the second pixel axis
00129 // from the CoordinateSystem, but leave the world axis intact.  This enables
00130 // us to still be able to make coordinate conversions for the first (roughly RA)
00131 // pixel axis.  Thus, CoordinateSystem::referenceValue would return a Vector of
00132 // length 3 (for RA, DEC and spectral), but CoordinateSystem::referencePixel
00133 // would return a vector length 2 (for RA and spectral).  
00134 // <p>
00135 // Now this CoordinateSystem has two Coordinates, a DirectionCoordinate and
00136 // a SpectralCoordinate, and let us state that that is the order in which
00137 // they exist in the CoordinateSystem (you can change them about if you wish);
00138 // they are coordinates number 0 and 1. The DirectionCoordinate has two axes
00139 // (RA and DEC) and the SpectralCoordinate has one axis. Only the
00140 // CoordinateSystem knows about removed axes, the DirectionCoordinate
00141 // itself is ignorant that it has been bisected. If you want to find
00142 // out what axis in the Coordinate system is where, you can use
00143 // the functions findPixelAxis or findWorldAxis.
00144 //
00145 // If we asked the former to find pixel axis 0, it would tell us that the
00146 // Coordinate number was 0 (the DirectionCoordinate) and that the axis in
00147 // that coordinate was 0 (the first axis in a DirectionCoordinate
00148 // is always longitude, the second always latitude).  If we asked it to find
00149 // pixel axis 1, it would tell us that the coordinate number was 1
00150 // (the SpectralCoordinate) and that the axis in that coordinate was 0
00151 // (there is only one axis in a SpectralCoordinate). If we asked for
00152 // pixelAxis 2 that would generate an error because our squashed image
00153 // only has 2 pixel axes. 
00154 //
00155 // Now, if we asked findWorldAxis similar questions,
00156 // it would tell us that worldAxis 0 in the CoordinateSystem can be found in
00157 // coordinate 0 (the DirectionCoordinate) in axis 0 of that DirectionCoordinate.
00158 // Similarly, worldAxis 1 in the CoordinateSystem (which has not been removed)
00159 // is in coordinate 0 (the DirectionCoordinate) in axis 1 of that 
00160 // Finally, worldAxis 2 in the CoordinateSystem is in coordinate 1 
00161 // (the SpectralCoordinate) in axis 0 of that SpectralCoordinate.
00162 // <p>
00163 // Other handy functions are pixelAxes and worldAxes.
00164 // These list the pixel and world axes in
00165 // the CoordinateSystem for the specified coordinate. Thus, if we asked
00166 // pixelAxes to find the pixel axes for coordinate 0 (the  DirectionCoordinate)
00167 // in the CoordinateSystem it would return a vector [0, -1]  indicating
00168 // the second axis of  the DirectionCoordinate has been removed.  However, 
00169 // the worldAxes function would return [0,1] as no world axis has been removed.
00170 // Similarly, if operated on coordinate 1 (the SpectralCoordinate), pixelAxes
00171 // would return [1] and worldAxes would return [2].
00172 //
00173 // Because you can transpose the CoordinateSystem about, you should NEVER ASSUME
00174 // ANYTHING except that the pixel axes of the CoordinateSystem map to the pixel
00175 // axes of the image when you first construct the image.
00176 //
00177 // <p>
00178 // SpectralCoordinate and DirectionCoordinate both have a (non-virtual) function
00179 // called <src>setReferenceConversion</src>.  This enables an extra conversion
00180 // layer so that conversion between pixel and world can go to a reference frame
00181 // other than the construction reference.    When you use the function
00182 // <src>convert</src>, these layers are active, but ONLY if the 
00183 // requested conversion is purely between pixel and world. For 
00184 // a SpectralCoordinate this must always be true (only has one axis)
00185 // but for the DirectionCoordinate  you might request a mixed
00186 // pixel/world conversion. In this case, the extra conversion layer
00187 // is ill-defined and not active (for the DirectionCoordinate part of it).
00188 // </synopsis>
00189 
00190 // <note role=caution>
00191 // All pixels coordinates are zero relative.
00192 // </note>
00193 
00194 // <example>
00195 // See the example in <linkto module=Coordinates>Coordinates.h</linkto>
00196 // and tCoordinateSystem.cc
00197 // </example>
00198 
00199 // <motivation>
00200 // Coordinate systems for images.
00201 // </motivation>
00202 //
00203 // <thrown>
00204 //   <li>  AipsError
00205 // </thrown>
00206 //
00207 // <todo asof="1997/01/13">
00208 //   <li> Undelete individual removed axes.
00209 //   <li> Non-integral pixel shifts/decimations in subimage operations?
00210 //   <li> Copy-on-write for efficiency?
00211 //   <li> Check if the classes are thread safe in general
00212 // </todo>
00213 //
00214 
00215 
00216 class CoordinateSystem : public Coordinate
00217 {
00218 public:
00219     // Default constructor.  This is an empty CoordinateSystem.
00220     CoordinateSystem();
00221 
00222     // Copying constructor (copy semantics)
00223     CoordinateSystem(const CoordinateSystem &other);
00224 
00225     // Assignment (copy semantics).
00226     CoordinateSystem &operator=(const CoordinateSystem &other);
00227 
00228     // Destructor
00229     virtual ~CoordinateSystem();
00230 
00231     // Add another Coordinate to this CoordinateSystem. This addition is done
00232     // by copying, so that if coord changes the change is NOT
00233     // reflected in the CoordinateSystem.
00234     void addCoordinate(const Coordinate &coord);
00235 
00236     // Transpose the CoordinateSystem so that world axis 0 is
00237     // newWorldOrder(0) and so on for all the other axes.
00238     // newPixelOrder works similarly. Normally you will give the
00239     // same transformation vector for both the world and pixel transformations,
00240     // however this is not required.
00241     void transpose(const Vector<Int> &newWorldOrder,
00242                    const Vector<Int> &newPixelOrder);
00243 
00244     // Find the world and pixel axis mappings to the supplied CoordinateSystem
00245     // from the current coordinate system. <src>False</src> is 
00246     // returned if either the supplied or current coordinate system, 
00247     // has no world axes (and a message recoverable with function
00248     // errorMessage indicating why).  Otherwise <src>True</src> is returned.
00249     // worldAxisMap(i) is the location of world axis <src>i</src> (from the
00250     // supplied CoordinateSystem, cSys, in the current CoordinateSystem.
00251     // worldAxisTranspose(i) is the location of world axis 
00252     // <src>i</src> (from the current CoordinateSystem) in the supplied 
00253     // CoordinateSystem, cSys.  The output vectors
00254     // are resized appropriately by this function.  A value of  -1 
00255     // in either vector means that the axis could not be found in the other
00256     // CoordinateSystem.  The vector <src>refChange</src> says
00257     // if the types are the same, is there a reference type change
00258     // (e.g. TOPO versus LSR for the SpectralCoordinate, 
00259     // or J2000 versus GALACTIC for DirectionCoordinate). Thus
00260     // if refChange(i) is True, it means world axis i in the
00261     // current CoordinateSystem was matched, but has a different
00262     // reference type to that of the supplied CoordinateSystem.
00263     // <group>
00264     Bool worldMap (Vector<Int>& worldAxisMap,
00265                    Vector<Int>& worldAxisTranspose,
00266                    Vector<Bool>& refChange,
00267                    const CoordinateSystem& cSys) const;
00268     Bool pixelMap (Vector<Int>& pixelAxisMap,
00269                    Vector<Int>& pixelAxisTranspose,
00270                    const CoordinateSystem& cSys) const;
00271     // </group>
00272 
00273     // Remove a world or pixel axis. When its value is required for forward or
00274     // backwards transformations, use <src>replacement</src>
00275     // <br>
00276     // When a world axis is removed, the corresponding pixel axis is removed
00277     // too, because it makes no sense having a pixel axis without world
00278     // coordinates.
00279     // <br>
00280     // Removing a pixel axis without removing the corresponding world axis
00281     // is, however, possible and meaningful. It can be used when e.g. a
00282     // frequency plane is taken from a cube. The plane has 2 pixel axes, but
00283     // the 3rd world axis can still describe the frequency coordinate.
00284     // See also the functions in  <linkto class=CoordinateUtil>CoordinateUtil</linkto>
00285     // for removing lists of pixel/world axes (tricky because they shift down)
00286     //
00287     // False is returned (an error in <src>errorMessage()</src> will be set)
00288     // if the axis is illegal, else returns True.
00289     // <group>
00290     Bool removeWorldAxis(uInt axis, Double replacement);
00291     Bool removePixelAxis(uInt axis, Double replacement);
00292     // </group>
00293 
00294     // Return a CoordinateSystem appropriate for a shift of origin
00295     // (the shift is subtracted from the reference pixel)
00296     // and change of increment (the increments are multipled
00297     // by the factor). Both vectors should be of length nPixelAxes(). 
00298     //
00299     // The newShape vector is only needed for the StokesCoordinate,
00300     // if any.  If this vector is of length zero, the new StokesCoordinate
00301     // is formed from all of the available input Stokes after application
00302     // of the shift and increment factor.    Otherwise,
00303     // the new Stokes axis length is equal to that specified after
00304     // appliction of the shift and increment and excess values 
00305     // discarded.    In addition, for any StokesCoordinate, the
00306     // shift and factor must be integer.  So <src>Int(value+0.5)</src>
00307     // is taken before they are used.
00308     // <group>
00309     CoordinateSystem subImage(const Vector<Float> &originShift,
00310                               const Vector<Float> &incrFac,
00311                               const Vector<Int>& newShape) const;
00312     void subImageInSitu (const Vector<Float> &originShift,
00313                          const Vector<Float> &incrFac,
00314                          const Vector<Int>& newShape);
00315     // </group>
00316 
00317     // Untranspose and undelete all axes. Does not undo the effects of
00318     // subimaging.
00319     void restoreOriginal();
00320 
00321     // Returns the number of Coordinates that this CoordinateSystem contains.
00322     // The order might be unrelated to the axis order through the results of
00323     // transposing and removing axes.
00324     uInt nCoordinates() const;
00325 
00326     // For a given Coordinate say where its world and pixel axes are in
00327     // this CoordinateSystem. The position in the returned Vector is its
00328     // axis number in the Coordinate, and its value is the axis
00329     // number in the CoordinateSystem. If the value is less than zero the axis
00330     // has been removed from this CoordinateSystem.
00331     //  <group>
00332     Vector<Int> worldAxes(uInt whichCoord) const;
00333     Vector<Int> pixelAxes(uInt whichCoord) const;
00334     // </group> 
00335 
00336     // Return the type of the given Coordinate.
00337     Coordinate::Type type(uInt whichCoordinate) const;
00338 
00339     // Returns the type of the given Coordinate as a string.
00340     String showType(uInt whichCoordinate) const;
00341 
00342     // Return the given Coordinate as a reference to the base
00343     // class object.
00344     const Coordinate& coordinate(uInt which) const;
00345 
00346     // Return the given Coordinate.
00347     // Throws an exception if retrieved as the wrong type.
00348     // The versions which take no parameters will return the
00349     // first (or in most cases only) coordinate of the requested type.
00350     // If no such coordinate exists, an exception is thrown.
00351     // <group>
00352     const LinearCoordinate    &linearCoordinate(uInt which) const;
00353     const DirectionCoordinate &directionCoordinate() const;
00354     const DirectionCoordinate &directionCoordinate(uInt which) const;
00355 
00356     const SpectralCoordinate &spectralCoordinate(uInt which) const;
00357     const SpectralCoordinate &spectralCoordinate() const;
00358     const StokesCoordinate  &stokesCoordinate() const;
00359 
00360     const StokesCoordinate  &stokesCoordinate(uInt which) const;
00361     const QualityCoordinate &qualityCoordinate(uInt which) const;
00362     const TabularCoordinate &tabularCoordinate(uInt which) const;
00363     // </group>
00364 
00365     // Replace one Coordinate with another. The mapping of the coordinate axes
00366     // to the CoordinateSystem axes is unchanged, therefore the number of world
00367     // and pixel axes must not be changed. You can, somewhat dangerously,
00368     // change the type of the coordinate however. For example, replace a 
00369     // SpectralCoordinate with a 1-D Linearcoordinate.  It is dangerous because
00370     // the world replacement values (see removeWorldAxis) have to be scaled.
00371     // The algorithm tries to find a scale factor between the old and new
00372     // units and applies it to the replacement values.  If it can't find
00373     // a scale factor (non-conformant units) then the reference value is
00374     // used for any world replacement values.  If the latter occurs,
00375     // it returns False, else True is returned.
00376     Bool replaceCoordinate(const Coordinate &newCoordinate, uInt whichCoordinate);
00377 
00378     // Find the Coordinate number that corresponds to the given type.
00379     // Since there might be more than one Coordinate of a given type you
00380     // can call this multiple times setting <src>afterCoord</src> to
00381     // the last value found. Returns -1 if a Coordinate of the desired
00382     // type is not found.
00383     Int findCoordinate(Coordinate::Type type, Int afterCoord = -1) const;
00384 
00385     // Given an axis number (pixel or world) in the CoordinateSystem,
00386     // find the corresponding coordinate number and axis in that Coordinate. 
00387     // The returned values are set to -1 if the axis does not exist.
00388     // <group>
00389     void findWorldAxis(Int &coordinate, Int &axisInCoordinate, 
00390                        uInt axisInCoordinateSystem) const;
00391     void findPixelAxis(Int &coordinate, Int &axisInCoordinate, 
00392                        uInt axisInCoordinateSystem) const;
00393     // </group>
00394 
00395     // Find the world axis for the given pixel axis in a CoordinateSystem.
00396     // Returns -1 if the world axis is unavailable (e.g. if it has been
00397     // removed).  
00398     Int pixelAxisToWorldAxis(uInt pixelAxis) const;
00399 
00400     // Find the pixel axis for the given world axis in a CoordinateSystem.
00401     // Returns -1 if the pixel axis is unavailable (e.g. if it has been
00402     // removed). 
00403     Int worldAxisToPixelAxis(uInt worldAxis) const;
00404 
00405     // Return the name of the record field in which the coordinate is stored.
00406     String coordRecordName(uInt which) const;
00407 
00408     // Returns <src>Coordinate::COORDSYS</src>
00409     virtual Coordinate::Type type() const;
00410 
00411     // Always returns "System"
00412     virtual String showType() const;
00413 
00414     // Sums the number of axes in the Coordinates that the CoordinateSystem
00415     // contains, allowing for removed axes.
00416     // <group>
00417     virtual uInt nPixelAxes() const;
00418     virtual uInt nWorldAxes() const;
00419     // </group>
00420 
00421 
00422     // Convert a pixel position to a world position or vice versa. Returns True
00423     // if the conversion succeeds, otherwise it returns <src>False</src> and
00424     // <src>errorMessage()</src> contains an error message. 
00425     // The input vector must be of length <src>nPixelAxes</src> or
00426     // <src>nWorldAxes</src>.  The output vector  is resized appropriately.
00427     // if <src>useConversionFrame</src>, if the coordinate has a conversion layer frame
00428     // (such as can be present in spectral and direction coordinates), it
00429     // is used. Else, the native frame is used for the conversion.
00430     // <group>
00431     virtual Bool toWorld(Vector<Double> &world, 
00432                          const Vector<Double> &pixel, Bool useConversionFrame=True) const;
00433     // This one throws an exception rather than returning False. After all, that's
00434     // what exceptions are for.
00435     virtual Vector<Double> toWorld(const Vector<Double> &pixel) const;
00436     virtual Bool toPixel(Vector<Double> &pixel, 
00437                          const Vector<Double> &world) const;
00438     // This one throws an exception rather than returning False.
00439     virtual Vector<Double> toPixel(const Vector<Double> &world) const;
00440     // </group>
00441 
00442     // convert a pixel "length" to a world "length"
00443     virtual Quantity toWorldLength(
00444         const Double nPixels,
00445         const uInt pixelAxis
00446     ) const;
00447 
00448     // This is provided as a convenience since it is a very commonly desired
00449     // operation through CoordinateSystem.  The output vector is resized.   
00450     Bool toWorld(Vector<Double> &world, const IPosition &pixel) const;
00451     Vector<Double> toWorld(const IPosition& pixel) const;
00452 
00453     // Batch up a lot of transformations. The first (most rapidly varying) axis
00454     // of the matrices contain the coordinates. Returns False if any conversion
00455     // failed  and  <src>errorMessage()</src> will hold a message.
00456     // The <src>failures</src> array (True for fail, False for success)
00457     // is the length of the number of conversions and
00458     // holds an error status for each conversion.  
00459     // <group>
00460     virtual Bool toWorldMany(Matrix<Double>& world,
00461                              const Matrix<Double>& pixel,
00462                              Vector<Bool>& failures) const;
00463     virtual Bool toPixelMany(Matrix<Double>& pixel,
00464                              const Matrix<Double>& world,
00465                              Vector<Bool>& failures) const;
00466     // </group>
00467 
00468 
00469     // Mixed pixel/world coordinate conversion.
00470     // <src>worldIn</src> and <src>worldAxes</src> are of length n<src>worldAxes</src>.
00471     // <src>pixelIn</src> and <src>pixelAxes</src> are of length nPixelAxes.
00472     // <src>worldAxes(i)=True</src> specifies you have given a world
00473     // value in <src>worldIn(i)</src> to convert to pixel.
00474     // <src>pixelAxes(i)=True</src> specifies you have given a pixel 
00475     // value in <src>pixelIn(i)</src> to convert to world.
00476     // You cannot specify the same axis via <src>worldAxes</src>
00477     // and pixelAxes.
00478     // Values in <src>pixelIn</src> are converted to world and
00479     // put into <src>worldOut</src> in the appropriate world axis
00480     // location.  Values in <src>worldIn</src> are copied to
00481     // <src>worldOut</src>.   
00482     // Values in <src>worldIn</src> are converted to pixel and
00483     // put into <src>pixelOut</src> in the appropriate pixel axis
00484     // location.  Values in <src>pixelIn</src> are copied to
00485     // <src>pixelOut</src>.  Vectors
00486     // <src>worldMin</src> and <src>worldMax</src> specify the range of the world
00487     // coordinate (in the world axis units of that world axis
00488     // in the coordinate system) being solved for in a mixed calculation 
00489     // for each world axis. They are only actually used for DirectionCoordinates
00490     // and for all other coordinates the relevant elements are ignored.
00491     // Functions <src>setWorldMixRanges, worldMixMin, worldMixMax</src> can be
00492     // used to compute and recover the world ranges.  If you don't know 
00493     // the values, use  functions <src>setDefaultWorldMixRanges, worldMixMin, worldMixMax</src>.
00494     // Removed axes are handled (for example, a removed pixel
00495     // axis with remaining corresponding world axis will
00496     // correctly be converted to world using the replacement
00497     // value).
00498     // Returns True if the conversion succeeds, otherwise it returns <src>False</src> and
00499     // <src>errorMessage()</src> contains an error message. The output vectors
00500     // are resized.
00501     virtual Bool toMix(Vector<Double>& worldOut,
00502                        Vector<Double>& pixelOut,
00503                        const Vector<Double>& worldIn,
00504                        const Vector<Double>& pixelIn,
00505                        const Vector<Bool>& worldAxes,                
00506                        const Vector<Bool>& pixelAxes,
00507                        const Vector<Double>& worldMin,
00508                        const Vector<Double>& worldMax) const; 
00509 
00510     // Compute and recover the world min and max ranges, for use in function <src>toMix</src>,
00511     // for  a lattice of the given shape (must be of length <src>nPixelAxes()</src>). 
00512     // Removed pixel axes (with remaining world axes are handled).  With
00513     // the retrieval functions, the output vectors are resized.  They return 
00514     // False if they fail (and then <src>setDefaultWorldMixRanges</src> generates the ranges)
00515     // with a reason in <src>errorMessage()</src>.
00516     // The <src>setDefaultWorldMixRanges</src> function
00517     // gives you  a useful default range if you don't know the shape.
00518     // The only Coordinate type for which these ranges are actually
00519     // used in <src>toMix</src> is DirectionCoordinate (because its coupled). For
00520     // the rest the functionality is provided but never used
00521     // by toMix.
00522     //<group>
00523     virtual Bool setWorldMixRanges (const IPosition& shape);
00524     virtual void setDefaultWorldMixRanges ();
00525     virtual Vector<Double> worldMixMin () const;
00526     virtual Vector<Double> worldMixMax () const;
00527     //</group>
00528 
00529     // Make absolute coordinates relative and vice-versa (relative
00530     // to the reference pixel/value).  The vectors must be of length
00531     // <src>nPixelAxes()</src> or <src>nWorldAxes()</src>
00532     //<group>
00533     virtual void makePixelRelative (Vector<Double>& pixel) const;
00534     virtual void makePixelAbsolute (Vector<Double>& pixel) const;
00535     virtual void makeWorldRelative (Vector<Double>& world) const;
00536     virtual void makeWorldAbsolute (Vector<Double>& world) const;
00537     //</group>
00538 
00539     // Make absolute coordinates relative and vice versa with respect
00540     // to the given reference value.  Add the other functions in this grouping
00541     // as needed.    The vectors must be of length
00542     // <src>nPixelAxes()</src> or <src>nWorldAxes()</src>
00543     //<group>
00544     virtual void makeWorldAbsoluteRef (Vector<Double>& world,
00545                                        const Vector<Double>& refVal) const;
00546     //</group>
00547 
00548     // Batch up a lot of absolute/relative transformations. 
00549     // Parameters as above  for 
00550     // <src>toWorldMany</src> and <src>toPixelMany</src>
00551     // <group>
00552     virtual void makePixelRelativeMany (Matrix<Double>& pixel) const;
00553     virtual void makePixelAbsoluteMany (Matrix<Double>& pixel) const;
00554     virtual void makeWorldRelativeMany (Matrix<Double>& world) const;
00555     virtual void makeWorldAbsoluteMany (Matrix<Double>& world) const;
00556     // </group>
00557 
00558 
00559     // General coordinate conversion.  Only works if no axes
00560     // have been removed and no axis reordering has occurred.
00561     // That is pixel axes and world axes are the same.
00562     //
00563     // Specify the input coordinate values, input units,
00564     // whether value is absolute (or relative). For output
00565     // specify units and abs/rel.  Units may be 'pix' and velocity consistent
00566     // units (e.g. m/s).  Specify doppler types if velocities
00567     // involved.   The pixel offsets allow for the input
00568     // and output pixel coordinates to be something other than 0-rel.  
00569     // If your pixel coordinates are 1-rel input and output, set the 
00570     // offsets to -1 and 1
00571     //
00572     // The Matrix interface lets you do many conversions efficiently.
00573     // Use <src>Matrix(nAxes, nConversions) </src> and 
00574     // <src>Matrix.column()=coordinate</src> or
00575     // <src>Matrix(axis, iConversion)</src> to get the order right.  
00576     //
00577     // These functions invoke <src>toMix</src>
00578     // so make sure you call <src>setWorldMixRanges</src>
00579     // first to set up the world ranges. 
00580     // <group>
00581     Bool convert (Vector<Double>& coordOut, 
00582                   const Vector<Double>& coordin,
00583                   const Vector<Bool>& absIn,
00584                   const Vector<String>& unitsIn,
00585                   MDoppler::Types dopplerIn,
00586                   const Vector<Bool>& absOut,
00587                   const Vector<String>& unitsOut,
00588                   MDoppler::Types dopplerOut,
00589                   Double pixInOffset = 0.0,
00590                   Double pixOutOffset = 0.0);
00591     Bool convert (Matrix<Double>& coordOut, 
00592                   const Matrix<Double>& coordIn,
00593                   const Vector<Bool>& absIn,
00594                   const Vector<String>& unitsIn,
00595                   MDoppler::Types dopplerIn,
00596                   const Vector<Bool>& absOut,
00597                   const Vector<String>& unitsOut,
00598                   MDoppler::Types dopplerOut,
00599                   Double pixInOffset = 0.0,
00600                   Double pixOutOffset = 0.0);
00601     // </group>
00602 
00603     // Return the requested attribute.
00604     // <group>
00605     virtual Vector<String> worldAxisNames() const;
00606     virtual Vector<Double> referencePixel() const;
00607     virtual Matrix<Double> linearTransform() const;
00608     virtual Vector<Double> increment() const;
00609     virtual Vector<Double> referenceValue() const;
00610     // </group>
00611 
00612     // Set the requested attribute.  Note that these just
00613     // change the internal values, they do not cause any recomputation.
00614     // <group>
00615     virtual Bool setWorldAxisNames(const Vector<String> &names);
00616     virtual Bool setReferencePixel(const Vector<Double> &refPix);
00617     virtual Bool setLinearTransform(const Matrix<Double> &xform);
00618     virtual Bool setIncrement(const Vector<Double> &inc);
00619     virtual Bool setReferenceValue(const Vector<Double> &refval);
00620     // </group>
00621 
00622     // Set/get the units. Adjust the increment and
00623     // reference value by the ratio of the old and new units. This implies that
00624     // the units must be known <linkto class=Unit>Unit</linkto> strings, and
00625     // that they must be compatible, e.g. they can't change from time to
00626     // length. If <src>throwException=True</src>, throw an exception rather than
00627     // returning False on failure.
00628     // <group>
00629     virtual Bool setWorldAxisUnits(const Vector<String> &units);
00630     Bool setWorldAxisUnits(const Vector<String> &units,
00631                            Bool throwException);
00632     virtual Vector<String> worldAxisUnits() const;
00633     // </group>
00634 
00635     // Comparison function. Any private Double data members are compared
00636     // with the specified fractional tolerance.  Don't compare on the specified 
00637     // pixel axes in the CoordinateSystem.  If the comparison returns
00638     // <src>False</src>, errorMessage() contains a message about why.
00639     // <group>
00640     virtual Bool near(const Coordinate& other, Double tol=1e-6) const;
00641     virtual Bool near(const Coordinate& other, 
00642                       const Vector<Int>& excludePixelAxes,
00643                       Double tol=1e-6) const;
00644     // </group>
00645 
00646     // This function compares this and the other coordinate system,
00647     // but ONLY for the non-removed pixel axes.   It is less strict
00648     // than near, which, for example, insists the number of coordinates
00649     // is the same in each CS
00650     Bool nearPixel (const CoordinateSystem& other, Double tol=1e-6) const;
00651 
00652 
00653     // Format a world value nicely through the
00654     // common format interface.  See <linkto class=Coordinate>Coordinate</linkto>
00655     // for basics.
00656     //
00657     // You specify a world value and its corresponding world axis in
00658     // the CoordinateSystem.   
00659     //
00660     // For the specified worldAxis, the coordinate
00661     // number in the CoordinateSystem is found and the actual derived Coordinate
00662     // class object for that number is created.  The arguments to the formatting 
00663     // function are then passed on to the formatter for that Coordinate. So
00664     // refer to the other derived Coordinate classes for specifics on the
00665     // formatting.
00666     virtual String format(
00667         String& units,
00668         Coordinate::formatType format,
00669         Double worldValue,
00670         uInt worldAxis,
00671         Bool isAbsolute=True,
00672         Bool showAsAbsolute=True,
00673         Int precision=-1, Bool usePrecForMixed=False
00674     ) const;
00675 
00676     // Miscellaneous information related to an observation, for example the
00677     // observation date.
00678     // <group>
00679     ObsInfo obsInfo() const;
00680     void setObsInfo(const ObsInfo &obsinfo);
00681     // </group>
00682 
00683     // Find the CoordinateSystem (you can safely caste the pointer to a CoordinateSystem)
00684     // for when we Fourier Transform ourselves.  This pointer 
00685     // must be deleted by the caller. Axes specifies which pixel axes of the Coordinate
00686     // System you wish to transform.   Shape specifies the shape of the image
00687     // associated with all the axes of the CoordinateSystem.  Currently you have
00688     // no control over the reference pixel, it is always shape/2.
00689     virtual Coordinate* makeFourierCoordinate (const Vector<Bool>& axes,
00690                                                const Vector<Int>& shape) const;
00691 
00692 
00693     // Save the CoordinateSystem into the supplied record using the supplied field name.
00694     // The field must not exist, otherwise <src>False</src> is returned.
00695     // If the CoordinateSystem is empty  <src>False</src> is also returned.
00696     // If <src>False</src> is returned, errorMessage() contains a message about why.   
00697     virtual Bool save(RecordInterface &container,
00698                     const String &fieldName) const;
00699 
00700     // Restore the CoordinateSystem from a record.  The <src>fieldName</src>
00701     // can be empty, in which case the CoordinateSystem is restored 
00702     // directly from the Record, rather than a subrecord of it.
00703     // A null pointer means that the restoration did not succeed - probably 
00704     // because fieldName doesn't exist or doesn't contain a CoordinateSystem.
00705     static CoordinateSystem *restore(const RecordInterface &container,
00706                                    const String &fieldName);
00707 
00708     // Make a copy of the CoordinateSystem using new. The caller is responsible for calling
00709     // delete.
00710     virtual Coordinate* clone() const;
00711 
00712     // Convert a CoordinateSystem to FITS, i.e. fill in ctype etc. In the record
00713     // the keywords are vectors, it is expected that the actual FITS code will
00714     // split them into scalars and upcase the names. Returns False if one of the
00715     // keywords is already taken.
00716     // 
00717     // If writeWCS is True, attempt to write the WCS convention (Greisen and
00718     // Calabretta "Representation of celestial coordinates in FITS"). 
00719     // Use <src>oneRelative=True</src> to convert zero-relative pixel coordinates to
00720     // one-relative FITS coordinates.
00721     //
00722     // prefix gives the prefix for the FITS keywords. E.g.,
00723     // if prefix="c" then crval, cdelt etc. 
00724     // if prefix="d" then drval, ddelt etc. 
00725     //# Much of the work in to/from fits should be moved to the individual
00726     //# classes.
00727     Bool toFITSHeader(RecordInterface &header, 
00728                       IPosition &shape,
00729                       Bool oneRelative, 
00730                       Char prefix = 'c', Bool writeWCS=True,
00731                       Bool preferVelocity=True, 
00732                       Bool opticalVelocity=True,
00733                       Bool preferWavelength=False,
00734                       Bool airWavelength=False) const;
00735 
00736     // Probably even if we return False we should set up the best linear
00737     // coordinate that we can.
00738     // Use oneRelative=True to convert one-relative FITS pixel coordinates to
00739     // zero-relative Casacore coordinates.
00740     // On output, <src>stokesFITSValue</src>
00741     // holds the FITS value of any unofficial Stokes (beam, optical depth,
00742     // spectral index) for the last unofficial value accessed (-1 if none).
00743     // The idea is that if the Stokes axis is of length one and holds an
00744     // unofficial value, you should drop the STokes axis and convert that
00745     // value to <src>ImageInfo::ImageTypes</src> with
00746     // <src>ImageInfo::imageTypeFromFITSValue</src>.
00747     // If on input, <src>stokesFITSValue</src> is positive, then a warning
00748     // is issued if any unofficial values are encountered.
00749     // Otherwise no warning is issued.
00750     //# cf comment in toFITS.
00751     static Bool fromFITSHeader(Int& stokesFITSValue, 
00752                                CoordinateSystem &coordsys, 
00753                                RecordInterface& recHeader,
00754                                const Vector<String>& header,
00755                                const IPosition& shape,
00756                                uInt which=0);
00757 
00758 // List all header information.  By default, the reference
00759 // values and pixel increments are converted to a "nice" unit before 
00760 // formatting (e.g. RA is  shown as HH:MM:SS.S).  
00761 // For spectral axes, both frequency and velocity information is listed. You
00762 // can specify what velocity definition you want with <src>velocityType</src>
00763 // If you wish, you can specify two shapes; a lattice and tile shape
00764 // (perhaps an image from which the CoordinateSystem came)
00765 // If you give (both of) these, they are included in the listing.  If you pass
00766 // in zero length <src>IPositions</src> then they are not included in
00767 // the listing.   If <src>postlocally=True</src> the formatted summary lines 
00768 // are written locally only to the sink, and then returned by the return value 
00769 // vector.
00770    Vector<String> list(LogIO& os, MDoppler::Types doppler,
00771                        const IPosition& latticeShape,
00772                        const IPosition& tileShape, Bool postLocally=False) const;
00773 
00774    // Does this coordinate system have a spectral axis?
00775    Bool hasSpectralAxis() const;
00776 
00777    // What number is the spectral axis?
00778    // If doWorld=True, the world axis number is returned.
00779    // Otherwise, the pixel axis number is returned.
00780    // Returns -1 if the spectral axis (world c.q. pixel) does not exist.
00781    Int spectralAxisNumber(Bool doWorld=False) const;
00782 
00783    // what number is the spectral coordinate?
00784     // Returns -1 if no spectral coordinate exists.
00785    Int spectralCoordinateNumber() const;
00786 
00787 
00788    // does this coordinate system have a polarizaion/stokes coordinate?
00789    Bool hasPolarizationCoordinate() const;
00790    Bool hasPolarizationAxis() const
00791      { return hasPolarizationCoordinate(); }
00792 
00793    // Given a stokes or polarization parameter, find the pixel location.
00794    // Note the client is responsible for any boundedness checks
00795    // (eg finite number of stokes in an image).
00796    Int stokesPixelNumber(const String& stokesString) const;
00797 
00798    // what is the number of the polarization/stokes coordinate?
00799    // Returns -1 if no stokes coordinate exists.
00800    Int polarizationCoordinateNumber() const;
00801 
00802    // What is the number of the polarization/stokes axis?
00803    // If doWorld=True, the world axis number is returned.
00804    // Otherwise, the pixel axis number is returned.
00805    // Returns -1 if the stokes axis (world c.q. pixel) does not exist.
00806    Int polarizationAxisNumber(Bool doWorld=False) const;
00807 
00808    // Does this coordinate system have a quality axis?
00809    Bool hasQualityAxis() const;
00810 
00811    // what number is the quality axis? Returns -1 if no quality axis exists.
00812    Int qualityAxisNumber() const;
00813 
00814    // what is the number of the quality coordinate?
00815    // Returns -1 if no quality coordinate exists.
00816    Int qualityCoordinateNumber() const;
00817 
00818    // Given a quality parameter, find the pixel location.
00819    // Note the client is responsible for any boundedness checks
00820    // (eg finite number of quality in an image).
00821    Int qualityPixelNumber(const String& qualityString) const;
00822 
00823    String qualityAtPixel(const uInt pixel) const;
00824 
00825    Int directionCoordinateNumber() const;
00826 
00827    Bool hasDirectionCoordinate() const;
00828 
00829    // Get the pixel axis numbers of the direction coordinate in this object.
00830    // The order of the returned axis numbers is always longitude axis first,
00831    // latitude axis second.
00832    Vector<Int> directionAxesNumbers() const;
00833 
00834    String stokesAtPixel(const uInt pixel) const;
00835 
00836    Int linearCoordinateNumber() const;
00837 
00838    Bool hasLinearCoordinate() const;
00839 
00840    Vector<Int> linearAxesNumbers() const;
00841 
00842    // Get the 0 based order of the minimal match strings specified in <src>order</src>.
00843    // If <src>requireAll</src> is True, checks are done to ensure that all axes in
00844    // the coordinate system are uniquely specified in <src>order</src>.
00845    // If <src>allowFriendlyNames</src> is True, the following (fully specified) strings
00846    // will match the specified axes:
00847    // "spectral" matches both "frequency" and "velocity".
00848    // "ra" matches "right ascension".
00849    Vector<Int> getWorldAxesOrder(Vector<String>& myNames, Bool requireAll,
00850                                  Bool allowFriendlyNames=False) const;
00851 
00852    // Is the abscissa in the DirectionCoordinate the longitude axis?
00853    // Throws exception if there is no DirectionCoordinate or if either of
00854    // the direction pixel axes have been removed.
00855    // For a normal direction coordinate, this will return True.
00856    Bool isDirectionAbscissaLongitude() const;
00857 
00858     // Set Spectral conversion layer of SpectralCoordinate in CoordinateSystem
00859     // so that pixel<->world go to the specified frequency system (a valid
00860     // MFrequency::Types string).  Returns False if frequency system invalid
00861     // or if no DirectionCoordinate or if cant get Date/Epoch.
00862     // <group>
00863     Bool setSpectralConversion (String& errorMsg, const String frequencySystem);
00864     // This version throws an exception rather than returning False.
00865     void setSpectralConversion (const String frequencySystem);
00866     //</group>
00867  
00868     // Set rest frequency of SpectralCoordinate in CoordinateSystem.
00869     // Unit must be consistent with Hz or m.
00870     // Returns False if invalid inputs (and CS not changed) and an error message.
00871     Bool setRestFrequency (String& errorMsg, const Quantity& freq);
00872 
00873 private:
00874     // Where we store copies of the coordinates we are created with.
00875     PtrBlock<Coordinate *> coordinates_p;
00876     
00877     // For coordinate[i] axis[j], 
00878     //    world_maps_p[i][j], if >=0 gives the location in the
00879     //                        input vector that maps to this coord/axis,
00880     //                        <0 means that the axis has been removed
00881     //    world_tmp_p[i] a temporary vector length coord[i]->nworldAxes()
00882     //    replacement_values_p[i][j] value to use for this axis if removed
00883     PtrBlock<Block<Int> *>     world_maps_p;
00884     PtrBlock<Vector<Double> *> world_tmps_p;
00885     PtrBlock<Vector<Double> *> world_replacement_values_p;
00886 
00887     // Same meanings as for the world*'s above.
00888     PtrBlock<Block<Int> *>     pixel_maps_p;
00889     PtrBlock<Vector<Double> *> pixel_tmps_p;
00890     PtrBlock<Vector<Double> *> pixel_replacement_values_p;
00891 
00892     // These temporaries all needed for the toMix function
00893     PtrBlock<Vector<Bool> *> worldAxes_tmps_p;
00894     PtrBlock<Vector<Bool> *> pixelAxes_tmps_p;
00895     PtrBlock<Vector<Double> *> worldOut_tmps_p;
00896     PtrBlock<Vector<Double> *> pixelOut_tmps_p;
00897     PtrBlock<Vector<Double> *> worldMin_tmps_p;
00898     PtrBlock<Vector<Double> *> worldMax_tmps_p;
00899 
00900     // Miscellaneous information about the observation associated with this
00901     // Coordinate System.
00902     ObsInfo obsinfo_p;
00903 
00904     const static String _class;
00905     static Mutex _mapInitMutex;
00906     static map<String, String> _friendlyAxisMap;
00907 
00908     static void _initFriendlyAxisMap();
00909 
00910     // Helper functions to group common code.
00911     Bool mapOne(Vector<Int>& worldAxisMap, 
00912                 Vector<Int>& worldAxisTranspose, 
00913                 Vector<Bool>& refChange,
00914                 const CoordinateSystem& cSys,
00915                 const CoordinateSystem& cSys2,
00916                 const uInt coord, const uInt coord2) const;
00917 
00918     void copy(const CoordinateSystem &other);
00919     void clear();
00920     Bool checkAxesInThisCoordinate(const Vector<Bool>& axes, uInt which) const;
00921 
00922    // Delete some pointer blocks
00923    void cleanUpSpecCoord (PtrBlock<SpectralCoordinate*>&  in,
00924                           PtrBlock<SpectralCoordinate*>&  out);
00925 
00926    // Delete temporary maps
00927    void deleteTemps (const uInt which);
00928 
00929     // Many abs/rel conversions
00930     // <group>
00931     void makeWorldAbsRelMany (Matrix<Double>& value, Bool toAbs) const;
00932     void makePixelAbsRelMany (Matrix<Double>& value, Bool toAbs) const;
00933     // </group>
00934 
00935     // Do subImage for Stokes
00936     StokesCoordinate stokesSubImage(const StokesCoordinate& sc, Int originShift, Int pixincFac,
00937                                     Int newShape) const;
00938 
00939     // Do subImage for Quality
00940     QualityCoordinate qualitySubImage(const QualityCoordinate& qc, Int originShift, Int pixincFac,
00941                                         Int newShape) const;
00942 
00943     // Strip out coordinates with all world and pixel axes removed
00944     CoordinateSystem stripRemovedAxes (const CoordinateSystem& cSys) const;
00945 
00946     //  All these functions are in support of the <src>list</src> function
00947     // <group>
00948     void listDirectionSystem(LogIO& os) const; 
00949     void listFrequencySystem(LogIO& os, MDoppler::Types velocityType) const;
00950     void listPointingCenter (LogIO& os) const;
00951     void getFieldWidths (LogIO& os, uInt& widthAxis, uInt& widthCoordType, 
00952                          uInt& widthCoordNumber, uInt& widthName,
00953                          uInt& widthProj, uInt& widthShape,
00954                          uInt& widthTile, uInt& widthRefValue,
00955                          uInt& widthRefPixel, uInt& widthInc,
00956                          uInt& widthUnits, Int& precRefValSci,
00957                          Int& precRefValFloat,  Int& precRefValRADEC,
00958                          Int& precRefPixFloat, Int& precIncSci, String& nameAxis,
00959                          String& nameCoordType, String& nameCoordNumber, String& nameName, String& nameProj,
00960                          String& nameShape, String& nameTile,
00961                          String& nameRefValue, String& nameRefPixel,
00962                          String& nameInc, String& nameUnits,
00963                          MDoppler::Types velocityType,
00964                          const IPosition& latticeShape, const IPosition& tileShape) const;
00965 
00966     void listHeader (LogIO& os, Coordinate* pc, uInt& widthAxis, uInt& widthCoordType,  uInt& widthCoordNumber,
00967                      uInt& widthName, uInt& widthProj,
00968                      uInt& widthShape, uInt& widthTile, uInt& widthRefValue,
00969                      uInt& widthRefPixel, uInt& widthInc, uInt& widthUnits,     
00970                      Bool findWidths, Int coordinate, Int axisInCoordinate, Int pixelAxis, 
00971                      Int precRefValSci, Int precRefValFloat, Int precRefValRADEC, Int precRefPixFloat,
00972                      Int precIncSci, const IPosition& latticeShape, const IPosition& tileShape) const;
00973     void listVelocity (LogIO& os,  Coordinate* pc, uInt widthAxis, 
00974                        uInt widthCoordType, uInt widthCoordNumber,
00975                        uInt& widthName, uInt widthProj,
00976                        uInt widthShape, uInt widthTile, uInt& widthRefValue,
00977                        uInt widthRefPixel, uInt& widthInc,  uInt& widthUnits,
00978                        Bool findWidths, Int axisInCoordinate, Int pixelAxis, 
00979                        MDoppler::Types velocityType, Int precRefValSci, Int precRefValFloat,
00980                        Int precRefValRADEC, Int precRefPixFloat, Int precIncSci) const;
00981     void clearFlags (LogIO& os) const;
00982     Bool velocityIncrement(Double& velocityInc,  SpectralCoordinate& sc,
00983                            MDoppler::Types velocityType, const String& velUnits) const;
00984     // </group>
00985 
00986     void _downcase(Vector<String>& vec) const
00987       { for (uInt i=0; i<vec.size(); ++i) vec[i].downcase(); }
00988 
00989 };
00990 
00991 } //# NAMESPACE CASACORE - END
00992 
00993 #endif
00994 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1