ImageRegrid.h

Go to the documentation of this file.
00001 //# ImageRegrid.h: Regrid Images
00002 //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
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 #ifndef IMAGES_IMAGEREGRID_H
00030 #define IMAGES_IMAGEREGRID_H
00031 
00032 #include <casacore/casa/aips.h>
00033 #include <casacore/casa/Arrays/Matrix.h>
00034 #include <casacore/casa/Arrays/Cube.h>
00035 #include <casacore/measures/Measures/MDirection.h>
00036 #include <casacore/measures/Measures/MFrequency.h>
00037 #include <casacore/scimath/Mathematics/Interpolate2D.h>
00038 #include <set>
00039 
00040 namespace casacore { //# NAMESPACE CASACORE - BEGIN
00041 
00042 template<class T> class MaskedLattice;
00043 template<class T> class ImageInterface;
00044 template<class T> class Lattice;
00045 template<class T> class LatticeIterator;
00046 template<class T> class Vector;
00047 
00048 class CoordinateSystem;
00049 class DirectionCoordinate;
00050 class Coordinate;
00051 class ObsInfo;
00052 class IPosition;
00053 class Unit;
00054 class ProgressMeter;
00055 
00056 // <summary>This regrids one image to match the coordinate system of another</summary>
00057 
00058 // <use visibility=export>
00059 
00060 // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
00061 // </reviewed>
00062 
00063 // <prerequisite>
00064 //   <li> <linkto class="ImageInterface">ImageInterface</linkto>
00065 //   <li> <linkto class="CoordinateSystem">CoordinateSystem</linkto>
00066 //   <li> <linkto class="Interpolate2D">Interpolate2D</linkto>
00067 //   <li> <linkto class="InterpolateArray1D">InterpolateArray1D</linkto>
00068 // </prerequisite>
00069 //
00070 // <etymology>
00071 //  Regrids, or resamples, images.  
00072 // </etymology>
00073 //
00074 // <synopsis>
00075 //  This class enables you to regrid one image to the coordinate
00076 //  system of another.    You can regrid any or all of the
00077 //  axes in the image.  A range of interpolation schemes are available.
00078 //
00079 //  It will cope with coordinate systems being in different orders
00080 //  (coordinate, world axes, pixel axes).  The basic approach is to
00081 //  make a mapping from the input to the output coordinate systems,
00082 //  but the output CoordinateSystem order is preserved in the output
00083 //  image.
00084 //
00085 //  Any DirectionCoordinate or LinearCoordinate holding exactly two axes
00086 //  is regridded in one pass with a 2-D interpolation scheme.
00087 //  All other axes are regridded in separate passes with a 1D interpolation 
00088 //  scheme.    This means that a LinearCoordinate holding say 3 axes
00089 //  where some of them are coupled will not be correctly regridded.
00090 //  StokesCoordinates cannot be  regridded.
00091 //
00092 //  Multiple passes are made through the data, and the output of 
00093 //  each pass is the input of the next pass.  The intermediate 
00094 //  images are stored as TempImages which may be in memory or 
00095 //  on disk, depending on their size.
00096 //
00097 //  It can also simply insert this image into that one via
00098 //  an integer shift.
00099 // </synopsis>
00100 //
00101 // <example>
00102 // 
00103 // <srcblock>
00104 // </srcblock>
00105 // </example>
00106 //
00107 // <motivation> 
00108 // A common image analysis need is to regrid images, e.g. to compare
00109 // images from different telescopes.
00110 // </motivation>
00111 //
00112 // <thrown>
00113 // <li> AipsError 
00114 // </thrown>
00115 //
00116 // <todo asof="1999/04/20">
00117 // </todo>
00118 
00119 template <class T> class ImageRegrid
00120 {
00121 public:
00122 
00123   // Default constructor
00124   ImageRegrid();
00125 
00126   // copy constructor (copy semantics)
00127   ImageRegrid(const ImageRegrid &other);
00128 
00129   // destructor
00130   ~ImageRegrid();
00131 
00132   // Assignment copy semantics)
00133   ImageRegrid<T>& operator=(const ImageRegrid& other);
00134 
00135   // Regrid inImage onto the grid specified by outImage.
00136   // If outImage has a writable mask, it will be updated in that 
00137   // output pixels at which the regridding failed will be masked bad (False)
00138   // and the pixel value set to zero. Otherwise the output mask is not changed.
00139   // Specify which pixel axes of outImage are to be
00140   // regridded.  The coordinate and axis order of outImage
00141   // is preserved, regardless of where the relevant coordinates
00142   // are in inImage.
00143   //
00144   // decimate only applies when replicate=False. it is
00145   // the coordinate grid computation decimation FACTOR
00146   // (e.g.  nCoordGrid ~ nIn / decimate). 0 means no decimation
00147   // (slowest and most accurate)
00148   void regrid(ImageInterface<T>& outImage, 
00149               typename Interpolate2D::Method method,
00150               const IPosition& whichOutPixelAxes,
00151               const ImageInterface<T>& inImage,
00152               Bool replicate=False, uInt decimate=0,
00153               Bool showProgress=False, Bool forceRegrid=False,
00154               Bool verbose=False);
00155 
00156 // Get and set the 2-D coordinate grid.  After a call to function <src>regrid</src>
00157 // in which coupled 2D coordinate (presently only DirectionCoordinate) is
00158 // regridded, this coordinate grid will be available.  It can be reused
00159 // via the <src>set2DCoordinateGrid</src> function for another like plane
00160 // (e.g. if you choose to regrid planes of a cube separately).   When you provide
00161 // the coordinate grid, it will no longer (for that 2D coordinate only) be
00162 // computed internally, which may save a lot of time.  Ordinarily, if you
00163 // regridded many planes of a cube in one call to regrid, the coordinate grid
00164 // is cached for you.   To trigger successive calls to regrid to go back to
00165 // internal computation, set zero length Cube and Matrix.  <src>gridMask</src>
00166 // is True for successfull coordinate conversions, and False otherwise.
00167 // <group>
00168   void get2DCoordinateGrid (Cube<Double>& grid, Matrix<Bool>& gridMask) const;
00169   void set2DCoordinateGrid (const Cube<Double>& grid, const Matrix<Bool>& gridMask, Bool notify=False);
00170 // </group>
00171 //
00172   // Inserts inImage into outImage.  The alignment is done by
00173   // placing the blc of inImage at the specified 
00174   // absolute pixel of the outImage (outPixelLocation).  If 
00175   // the outPixelLocation vector is of zero length, then the images 
00176   // are aligned by their reference pixels.  Only integral shifts are done
00177   // in the aligment process. If outImage has a mask,  it will be updated.
00178   // Returns False if no overlap of images, in which case the
00179   // output is not updated.
00180   Bool insert(ImageInterface<T>& outImage,
00181               const Vector<Double>& outPixelLocation,
00182               const ImageInterface<T>& inImage);
00183 
00184   // Print out useful debugging information (level 0 is none,
00185   // 1 is some, 2 is too much)
00186   void showDebugInfo(Int level=0) {itsShowLevel = level;};
00187 
00188   // Enable/disable Measures Reference conversions
00189   void disableReferenceConversions(Bool disable=True) {itsDisableConversions = disable;};
00190 
00191   // Helper function.  We are regridding from cSysFrom to cSysTo for the
00192   // specified pixel axes of cSyFrom. This function returns a CoordinateSystem which,
00193   // for the pixel axes being regridded, copies the coordinates from cSysTo
00194   // (if coordinate type present in cSysTo) or cSysFrom (coordinate
00195   // type not present in cSysTo).
00196   // For the axes not being regridded, it copies the coordinates from
00197   // cSysFrom.  This helps you build the cSys for function regrid.
00198   // The ObsInfo from cSysFrom is copied to the output CoordinateSystem.
00199   // If inShape has one or more elements it represenents the size of the
00200   // image to be regridded. It this must have the same number of elements
00201   // as the number of pixel axes in <src>cSysFrom</src>. If any of the values
00202   // are unity (ie the axes are degenerate), and the corresponding axis in <src>csysFrom</src> is the only
00203   // axis in its corresponding coordinate, this coordinate will not be replaced
00204   // even if the axis is specified in <src>axes</src>.
00205   // Upon return, <src>coordsToBeRegridded</src> will contain a list of the coordinates that will
00206   // be regridded.
00207   static CoordinateSystem makeCoordinateSystem(
00208                   LogIO& os,
00209                   std::set<Coordinate::Type>& coordsToBeRegridded,
00210                   const CoordinateSystem& cSysTo,
00211                   const CoordinateSystem& cSysFrom,
00212                   const IPosition& axes,
00213                   const IPosition& inShape=IPosition(),
00214                   Bool giveStokesWarning=True
00215   );
00216 
00217  private:
00218 
00219   Int itsShowLevel;
00220   Bool itsDisableConversions;
00221 //
00222   Cube<Double> its2DCoordinateGrid;
00223   Matrix<Bool> its2DCoordinateGridMask;
00224 //
00225   Cube<Double> itsUser2DCoordinateGrid;
00226   Matrix<Bool> itsUser2DCoordinateGridMask;
00227   Bool itsNotify;
00228 //  
00229   // Check shape and axes.  Exception if no good.  If pixelAxes
00230   // of length 0, set to all axes according to shape
00231   void _checkAxes(IPosition& outPixelAxes,
00232                   const IPosition& inShape,
00233                   const IPosition& outShape,
00234                   const Vector<Int>& pixelAxisMap,
00235                   const CoordinateSystem& outCoords,
00236                   Bool verbose);
00237 
00238   // Find maps between coordinate systems
00239   void findMaps (uInt nDim, 
00240                  Vector<Int>& pixelAxisMap1,
00241                  Vector<Int>& pixelAxisMap2,
00242                  const CoordinateSystem& inCoords,
00243                  const CoordinateSystem& outCoords) const;
00244 
00245   // Find scale factor to conserve flux 
00246    Double findScaleFactor(const Unit& units, 
00247                           const CoordinateSystem& inCoords, 
00248                           const CoordinateSystem& outCoords, 
00249                           Int inCoordinate, Int outCoordinate,
00250                           LogIO& os) const;
00251 
00252   // Regrid one Coordinate
00253    void _regridOneCoordinate (LogIO& os, IPosition& outShape2,
00254                               Vector<Bool>& doneOutPixelAxes,
00255                               MaskedLattice<T>* &finalOutPtr,  
00256                               MaskedLattice<T>* &inPtr,   
00257                               MaskedLattice<T>* &outPtr,  
00258                               CoordinateSystem& outCoords,
00259                               const CoordinateSystem& inCoords,
00260                               Int outPixelAxis,
00261                               const ImageInterface<T>& inImage,
00262                               const IPosition& outShape,
00263                               Bool replicate, uInt decimate,
00264                               Bool outIsMasked, Bool showProgress,
00265                               Bool forceRegrid, 
00266                               typename Interpolate2D::Method method,
00267                               Bool verbose);
00268 
00269   // Regrid  DirectionCoordinate or 2-axis LinearCoordinate
00270    void regridTwoAxisCoordinate  (LogIO& os, MaskedLattice<T>& outLattice,
00271                          const MaskedLattice<T>& inLattice,
00272                          const Unit& imageUnit, 
00273                          const CoordinateSystem& inCoords,
00274                          const CoordinateSystem& outCoords,
00275                          Int inCoordinate, Int outCoordinate,
00276                          const Vector<Int> inPixelAxes,
00277                          const Vector<Int> outPixelAxes,
00278                          const Vector<Int> pixelAxisMap1,  
00279                          const Vector<Int> pixelAxisMap2,
00280                          typename Interpolate2D::Method method,
00281                          Bool replicate, uInt decimate,
00282                          Bool showProgress);
00283 
00284   // Make regridding coordinate grid for this cursor.
00285   void make2DCoordinateGrid (LogIO& os, Bool& allFail, Bool&missedIt,
00286                              Double& minInX, Double& minInY, 
00287                              Double& maxInX, Double& maxInY,
00288                              Cube<Double>& in2DPos,
00289                              Matrix<Bool>& succeed,
00290                              const CoordinateSystem& inCoords,
00291                              const CoordinateSystem& outCoords,
00292                              Int inCoordinate, Int outCoordinate,
00293                              uInt xInAxis, uInt yInAxis,
00294                              uInt xOutAxis, uInt yOutAxis,
00295                              const IPosition& inPixelAxes,
00296                              const IPosition& outPixelAxes,
00297                              const IPosition& inShape,
00298                              const IPosition& outPos,
00299                              const IPosition& cursorShape,
00300                              uInt decimate=0);
00301 
00302   // Make replication coordinate grid for this cursor
00303    void make2DCoordinateGrid (Cube<Double>& in2DPos,
00304                               Double& minInX, Double& minInY, 
00305                               Double& maxInX, Double& maxInY,
00306                               const Vector<Double>& pixelScale, 
00307                               uInt xInAxis, uInt yInAxis,
00308                               uInt xOutAxis, uInt yOutAxis,
00309                               uInt xInCorrAxis, uInt yInCorrAxis,
00310                               uInt xOutCorrAxis, uInt yOutCorrAxis,
00311                               const IPosition& outPos, const IPosition& cursorShape);
00312 
00313   // Make regridding coordinate grid for this axis
00314    void make1DCoordinateGrid (Block<Float>& xOut,
00315                               Vector<Bool>& failed,
00316                               Bool& allFailed,
00317                               Bool& allGood,
00318                               const Coordinate& inCoord,
00319                               const Coordinate& outCoord,
00320                               Int inAxisInCoordinate,
00321                               Int outAxisInCoordinate,
00322                               MFrequency::Convert& machine,
00323                               Bool useMachine);
00324 
00325 
00326   // Make replication coordinate grid for this axis
00327    void make1DCoordinateGrid (Block<Float>& xOut, Float pixelScale) const;
00328 
00329   // Regrid 1 axis
00330   void regrid1D (MaskedLattice<T>& outLattice,
00331                  const MaskedLattice<T>& inLattice,
00332                  const Coordinate& inCoord,
00333                  const Coordinate& outCoord,
00334                  const Vector<Int>& inPixelAxes,
00335                  const Vector<Int>& outPixelAxes,
00336                  Int inAxisInCoordinate,
00337                  Int outAxisInCoordinate,
00338                  const Vector<Int> pixelAxisMap,
00339                  typename Interpolate2D::Method method,
00340                  MFrequency::Convert& machine,
00341                  Bool replicate,
00342                  Bool useMachine, Bool showProgress);
00343 
00344 //
00345    void regrid2DMatrix(Lattice<T>& outCursor,
00346                        LatticeIterator<Bool>*& outMaskIterPtr,
00347                        const Interpolate2D& interp,  
00348                                     ProgressMeter*& pProgress,
00349                                     Double& iPix,
00350                                     uInt nDim,
00351                                     uInt xInAxis, uInt yInAxis,
00352                                     uInt xOutAxis, uInt yOutAxis,
00353                                     Double scale,
00354                                     Bool inIsMasked, Bool outIsMasked,
00355                                     const IPosition& outPos,
00356                                     const IPosition& outCursorShape,
00357                                     const IPosition& inChunkShape,
00358                                     const IPosition& inChunkBlc,
00359                                     const IPosition& pixelAxisMap2,
00360                                     Array<T>& inDataChunk,
00361                                     Array<Bool>*& inMaskChunkPtr,
00362                                     const Cube<Double>& pix2DPos,
00363                                     const Matrix<Bool>& succeed);
00364 
00365    void findXYExtent (Bool& missedIt, Bool& allFailed,
00366                       Double& minInX, Double& minInY,
00367                       Double& maxInX, Double& maxInY,
00368                       Cube<Double>& in2DPos,
00369                       const Matrix<Bool>& succeed,
00370                       uInt xInAxis, uInt yInAxis,
00371                       uInt xOutAxis, uInt yOutAxis,
00372                       const IPosition& outPos,
00373                       const IPosition& outCursorShape,
00374                       const IPosition& inShape);
00375 //
00376    Bool minmax(Double &minX, Double &maxX, Double& minY, Double& maxY,
00377                const Array<Double> &xData,
00378                const Array<Double> &yData,
00379                const Array<Bool>& mask);
00380 };
00381 
00382 //# Declare extern templates for often used types.
00383 #ifdef AIPS_CXX11
00384   extern template class ImageRegrid<Float>;
00385 #endif
00386 
00387 } //# NAMESPACE CASACORE - END
00388 
00389 #ifndef CASACORE_NO_AUTO_TEMPLATES
00390 #include <casacore/images/Images/ImageRegrid.tcc>
00391 #endif //# CASACORE_NO_AUTO_TEMPLATES
00392 #endif
00393 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1