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