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

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1