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