WProjectFT.h

Go to the documentation of this file.
00001 //# WProjectFT.h: Definition for WProjectFT
00002 //# Copyright (C) 2003-2016
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 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 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 adressed 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 SYNTHESIS_TRANSFORM2_WPROJECTFT_H
00030 #define SYNTHESIS_TRANSFORM2_WPROJECTFT_H
00031 
00032 #include <synthesis/TransformMachines2/FTMachine.h>
00033 #include <casa/Arrays/Matrix.h>
00034 #include <scimath/Mathematics/FFTServer.h>
00035 #include <msvis/MSVis/VisBuffer2.h>
00036 #include <images/Images/ImageInterface.h>
00037 #include <images/Images/ImageInterface.h>
00038 #include <casa/Containers/Block.h>
00039 #include <casa/Arrays/Array.h>
00040 #include <casa/Arrays/Vector.h>
00041 #include <casa/Arrays/Matrix.h>
00042 #include <scimath/Mathematics/ConvolveGridder.h>
00043 #include <lattices/Lattices/LatticeCache.h>
00044 #include <lattices/Lattices/ArrayLattice.h>
00045 #include <ms/MeasurementSets/MSColumns.h>
00046 #include <measures/Measures/Measure.h>
00047 #include <measures/Measures/MDirection.h>
00048 #include <measures/Measures/MPosition.h>
00049 #include <coordinates/Coordinates/DirectionCoordinate.h>
00050 
00051 namespace casa { //# NAMESPACE CASA - BEGIN
00052 
00053 
00054 template <class K, class V> class SimpleOrderedMap;
00055 template <class T> class PtrBlock;
00056 template <class T> class CountedPtr;
00057 namespace refim { //# namespace for imaging refactor
00058 class   WPConvFunc; 
00059 
00060 // <summary>  An FTMachine for Gridded Fourier transforms </summary>
00061 
00062 // <use visibility=export>
00063 
00064 // <reviewed reviewer="" date="" tests="" demos="">
00065 
00066 // <prerequisite>
00067 //   <li> <linkto class=FTMachine>FTMachine</linkto> module
00068 //   <li> <linkto class=SkyEquation>SkyEquation</linkto> module
00069 //   <li> <linkto class=VisBuffer>VisBuffer</linkto> module
00070 // </prerequisite>
00071 //
00072 // <etymology>
00073 // FTMachine is a Machine for Fourier Transforms. WProjectFT does
00074 // Grid-based Fourier transforms.
00075 // </etymology>
00076 //
00077 // <synopsis> 
00078 // The <linkto class=SkyEquation>SkyEquation</linkto> needs to be able
00079 // to perform Fourier transforms on visibility data. WProjectFT
00080 // allows efficient Fourier Transform processing using a 
00081 // <linkto class=VisBuffer>VisBuffer</linkto> which encapsulates
00082 // a chunk of visibility (typically all baselines for one time)
00083 // together with all the information needed for processing
00084 // (e.g. UVW coordinates).
00085 //
00086 // Gridding and degridding in WProjectFT are performed using a
00087 // novel sort-less algorithm. In this approach, the gridded plane is
00088 // divided into small patches, a cache of which is maintained in memory
00089 // using a general-purpose <linkto class=LatticeCache>LatticeCache</linkto> class. As the (time-sorted)
00090 // visibility data move around slowly in the Fourier plane, patches are
00091 // swapped in and out as necessary. Thus, optimally, one would keep at
00092 // least one patch per baseline.  
00093 //
00094 // A grid cache is defined on construction. If the gridded uv plane is smaller
00095 // than this, it is kept entirely in memory and all gridding and
00096 // degridding is done entirely in memory. Otherwise a cache of tiles is
00097 // kept an paged in and out as necessary. Optimally the cache should be
00098 // big enough to hold all polarizations and frequencies for all
00099 // baselines. The paging rate will then be small. As the cache size is
00100 // reduced below this critical value, paging increases. The algorithm will
00101 // work for only one patch but it will be very slow!
00102 //
00103 // This scheme works well for arrays having a moderate number of
00104 // antennas since the saving in space goes as the ratio of
00105 // baselines to image size. For the ATCA, VLBA and WSRT, this ratio is
00106 // quite favorable. For the VLA, one requires images of greater than
00107 // about 200 pixels on a side to make it worthwhile.
00108 //
00109 // The FFT step is done plane by plane for images having less than
00110 // 1024 * 1024 pixels on each plane, and line by line otherwise.
00111 //
00112 // The gridding and degridding steps are implemented in Fortran
00113 // for speed. In gridding, the visibilities are added onto the
00114 // grid points in the neighborhood using a weighting function.
00115 // In degridding, the value is derived by a weight summ of the
00116 // same points, using the same weighting function.
00117 // </synopsis> 
00118 //
00119 // <example>
00120 // See the example for <linkto class=SkyModel>SkyModel</linkto>.
00121 // </example>
00122 //
00123 // <motivation>
00124 // Define an interface to allow efficient processing of chunks of 
00125 // visibility data
00126 // </motivation>
00127 //
00128 // <todo asof="97/10/01">
00129 // <ul> Deal with large VLA spectral line case 
00130 // </todo>
00131 
00132 class WProjectFT : public FTMachine {
00133 public:
00134 
00135   // Constructor: cachesize is the size of the cache in words
00136   // (e.g. a few million is a good number), tilesize is the
00137   // size of the tile used in gridding (cannot be less than
00138   // 12, 16 works in most cases). 
00139   // <group>
00140   WProjectFT(
00141            Int nFacets, Long cachesize, Int tilesize=16, 
00142            Bool usezero=True, Bool useDoublePrec=False, const Double minW=-1.0, const Double maxW=-1.0, const Double rmsW=-1.0);
00143   //Constructor without tangent direction
00144   WProjectFT(Int nFacets, MPosition mLocation,
00145              Long cachesize, Int tilesize=16, 
00146              Bool usezero=True, Float padding=1.0, Bool useDoublePrec=False, const Double minW=-1.0, const Double maxW=-1.0, const Double rmsW=-1.0);
00147   //Deprecated no longer need ms in constructor
00148   WProjectFT(
00149              Int nFacets, MDirection mTangent, MPosition mLocation,
00150              Long cachesize, Int tilesize=16, 
00151              Bool usezero=True, Float padding=1.0, Bool useDoublePrec=False, const Double minW=-1.0, const Double maxW=-1.0, const Double rmsW=-1.0);
00152   // </group>
00153 
00154   // Construct from a Record containing the WProjectFT state
00155   WProjectFT(const RecordInterface& stateRec);
00156 
00157   // Copy constructor
00158   WProjectFT(const WProjectFT &other);
00159 
00160   // Assignment operator
00161   WProjectFT &operator=(const WProjectFT &other);
00162 
00163   ~WProjectFT();
00164 
00165   //clone to FTMachine pointer
00166   virtual FTMachine* cloneFTM();
00167   // Initialize transform to Visibility plane using the image
00168   // as a template. The image is loaded and Fourier transformed.
00169   void initializeToVis(ImageInterface<Complex>& image,
00170                        const vi::VisBuffer2& vb);
00171   
00172   // Finalize transform to Visibility plane: flushes the image
00173   // cache and shows statistics if it is being used.
00174   void finalizeToVis();
00175 
00176   // Initialize transform to Sky plane: initializes the image
00177   void initializeToSky(ImageInterface<Complex>& image,  Matrix<Float>& weight,
00178                        const vi::VisBuffer2& vb);
00179   // Finalize transform to Sky plane: flushes the image
00180   // cache and shows statistics if it is being used. DOES NOT
00181   // DO THE FINAL TRANSFORM!
00182   void finalizeToSky();
00183 
00184   // Get actual coherence from grid by degridding
00185   void get(vi::VisBuffer2& vb, Int row=-1);
00186 
00187   // Put coherence to grid by gridding.
00188   void put(const vi::VisBuffer2& vb, Int row=-1, Bool dopsf=False,
00189            FTMachine::Type type=FTMachine::OBSERVED);
00190 
00191   // Make the entire image
00192   void makeImage(FTMachine::Type type,
00193                  vi::VisibilityIterator2& vi,
00194                  ImageInterface<Complex>& image,
00195                  Matrix<Float>& weight);
00196   
00197   // Get the final image: do the Fourier transform and
00198   // grid-correct, then optionally normalize by the summed weights
00199   ImageInterface<Complex>& getImage(Matrix<Float>&, Bool normalize=True);
00200   virtual void normalizeImage(Lattice<Complex>& /*skyImage*/,
00201                               const Matrix<Double>& /*sumOfWts*/,
00202                               Lattice<Float>& /*sensitivityImage*/,
00203                               Bool /*fftNorm*/)
00204     {throw(AipsError("WProjectFT::normalizeImage() called"));}
00205  
00206   // Get the final weights image
00207   void getWeightImage(ImageInterface<Float>&, Matrix<Float>&);
00208 
00209   // Save and restore the WProjectFT to and from a record
00210   Bool toRecord(String& error, RecordInterface& outRec, 
00211                 Bool withImage=False, const String diskimage="");
00212   Bool fromRecord(String& error, const RecordInterface& inRec);
00213   
00214   // Can this FTMachine be represented by Fourier convolutions?
00215   Bool isFourier() {return True;}
00216 
00217 
00218   // Return name of this machine
00219 
00220   String name() const;
00221 
00222   // Copy convolution function etc to another FT machine
00223   // necessary if ft and ift are distinct but can share convfunctions
00224 
00225   void setConvFunc(CountedPtr<WPConvFunc>& pbconvFunc);
00226   CountedPtr<WPConvFunc>& getConvFunc();
00227   virtual void setMiscInfo(const Int qualifier){(void)qualifier;};
00228   virtual void ComputeResiduals(vi::VisBuffer2& /*vb*/, Bool /*useCorrected*/) {};
00229   //Helper function to calculate min, max, rms of W in the data set
00230   static void wStat(vi::VisibilityIterator2& vi, Double& minW, Double& maxW, Double& rmsW); 
00231 
00232 protected:
00233 
00234   // Padding in FFT
00235   Float padding_p;
00236 
00237   Int nint(Double val) {return Int(floor(val+0.5));};
00238 
00239   // Find the convolution function
00240   void findConvFunction(const ImageInterface<Complex>& image,
00241                         const vi::VisBuffer2& vb);
00242 
00243   Int nWPlanes_p;
00244 
00245   // Get the appropriate data pointer
00246   Array<Complex>* getDataPointer(const IPosition&, Bool);
00247 
00248   void ok();
00249 
00250   void init();
00251 
00252   void prepGridForDegrid();
00253   // Is this record on Grid? check both ends. This assumes that the
00254   // ends bracket the middle
00255   //Bool recordOnGrid(const vi::VisBuffer2& vb, Int rownr) const;
00256 
00258   void   findGridSector(const Int& nxp, const Int& nyp, const Int& ixsub, const Int& iysub, const Int& minx, const Int& miny, const Int& icounter, Int& x0, Int& y0, Int& nxsub, Int& nysub, const Bool linear); 
00259 
00260 
00261   // Image cache
00262   LatticeCache<Complex> * imageCache;
00263 
00264   // Sizes
00265   Long cachesize;
00266   Int tilesize;
00267 
00268   // Gridder
00269   ConvolveGridder<Double, Complex>* gridder;
00270 
00271   // Is this tiled?
00272   Bool isTiled;
00273 
00274   // Array lattice
00275   CountedPtr<Lattice<Complex> > arrayLattice;
00276 
00277   // Lattice. For non-tiled gridding, this will point to arrayLattice,
00278   //  whereas for tiled gridding, this points to the image
00279   CountedPtr<Lattice<Complex> > lattice;
00280 
00281   Float maxAbsData;
00282 
00283   // Useful IPositions
00284   IPosition centerLoc, offsetLoc;
00285 
00286   // Image Scaling and offset
00287   Vector<Double> uvScale, uvOffset;
00288   Double savedWScale_p;
00289 
00290 
00291   // Grid/degrid zero spacing points?
00292   Bool usezero_p;
00293 
00294   Cube<Complex> convFunc;
00295   Int convSampling;
00296   Int convSize;
00297   Vector<Int> convSupport;
00298 
00299   Vector<Int> convSizes_p;
00300 
00301 
00302   Int wConvSize;
00303 
00304   Int lastIndex_p;
00305 
00306   Int getIndex(const ROMSPointingColumns& mspc, const Double& time,
00307                const Double& interval);
00308 
00309   String machineName_p;
00310 
00311   CountedPtr<WPConvFunc> wpConvFunc_p;
00312   Double timemass_p, timegrid_p, timedegrid_p;
00313   Double minW_p, maxW_p, rmsW_p;
00314 };
00315 } // end namespace refim
00316 } //# NAMESPACE CASA - END
00317 
00318 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1