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