00001 //# GridFT.h: Definition for GridFT 00002 //# Copyright (C) 1996-2012 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 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_GRIDFT_H 00030 #define SYNTHESIS_TRANSFORM2_GRIDFT_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 00046 00047 namespace casa { //# NAMESPACE CASA - BEGIN 00048 00049 class UVWMachine; 00050 namespace vi { class VisBuffer2;} 00051 namespace refim { //#namespace for imaging refactor 00052 // <summary> An FTMachine for Gridded Fourier transforms </summary> 00053 00054 // <use visibility=export> 00055 00056 // <reviewed reviewer="" date="" tests="" demos=""> 00057 00058 // <prerequisite> 00059 // <li> <linkto class=FTMachine>FTMachine</linkto> module 00060 // <li> <linkto class=SkyEquation>SkyEquation</linkto> module 00061 // <li> <linkto class=VisBuffer>VisBuffer</linkto> module 00062 // </prerequisite> 00063 // 00064 // <etymology> 00065 // FTMachine is a Machine for Fourier Transforms. GridFT does 00066 // Grid-based Fourier transforms. 00067 // </etymology> 00068 // 00069 // <synopsis> 00070 // The <linkto class=SkyEquation>SkyEquation</linkto> needs to be able 00071 // to perform Fourier transforms on visibility data. GridFT 00072 // allows efficient Fourier Transform processing using a 00073 // <linkto class=VisBuffer>VisBuffer</linkto> which encapsulates 00074 // a chunk of visibility (typically all baselines for one time) 00075 // together with all the information needed for processing 00076 // (e.g. UVW coordinates). 00077 // 00078 // Gridding and degridding in GridFT are performed using a 00079 // novel sort-less algorithm. In this approach, the gridded plane is 00080 // divided into small patches, a cache of which is maintained in memory 00081 // using a general-purpose <linkto class=LatticeCache>LatticeCache</linkto> class. As the (time-sorted) 00082 // visibility data move around slowly in the Fourier plane, patches are 00083 // swapped in and out as necessary. Thus, optimally, one would keep at 00084 // least one patch per baseline. 00085 // 00086 // A grid cache is defined on construction. If the gridded uv plane is smaller 00087 // than this, it is kept entirely in memory and all gridding and 00088 // degridding is done entirely in memory. Otherwise a cache of tiles is 00089 // kept an paged in and out as necessary. Optimally the cache should be 00090 // big enough to hold all polarizations and frequencies for all 00091 // baselines. The paging rate will then be small. As the cache size is 00092 // reduced below this critical value, paging increases. The algorithm will 00093 // work for only one patch but it will be very slow! 00094 // 00095 // This scheme works well for arrays having a moderate number of 00096 // antennas since the saving in space goes as the ratio of 00097 // baselines to image size. For the ATCA, VLBA and WSRT, this ratio is 00098 // quite favorable. For the VLA, one requires images of greater than 00099 // about 200 pixels on a side to make it worthwhile. 00100 // 00101 // The FFT step is done plane by plane for images having less than 00102 // 1024 * 1024 pixels on each plane, and line by line otherwise. 00103 // 00104 // The gridding and degridding steps are implemented in Fortran 00105 // for speed. In gridding, the visibilities are added onto the 00106 // grid points in the neighborhood using a weighting function. 00107 // In degridding, the value is derived by a weight summ of the 00108 // same points, using the same weighting function. 00109 // </synopsis> 00110 // 00111 // <example> 00112 // See the example for <linkto class=SkyModel>SkyModel</linkto>. 00113 // </example> 00114 // 00115 // <motivation> 00116 // Define an interface to allow efficient processing of chunks of 00117 // visibility data 00118 // </motivation> 00119 // 00120 // <todo asof="97/10/01"> 00121 // <ul> Deal with large VLA spectral line case 00122 // </todo> 00123 00124 class GridFT : public FTMachine { 00125 public: 00126 00127 // Constructor: cachesize is the size of the cache in words 00128 // (e.g. a few million is a good number), tilesize is the 00129 // size of the tile used in gridding (cannot be less than 00130 // 12, 16 works in most cases), and convType is the type of 00131 // gridding used (SF is prolate spheriodal wavefunction, 00132 // and BOX is plain box-car summation). mLocation is 00133 // the position to be used in some phase rotations. If 00134 // mTangent is specified then the uvw rotation is done for 00135 // that location iso the image center. 00136 // <group> 00137 GridFT(); 00138 GridFT(Long cachesize, Int tilesize, String convType="SF", 00139 Float padding=1.0, Bool usezero=True, Bool useDoublePrec=False); 00140 GridFT(Long cachesize, Int tilesize, String convType, 00141 MPosition mLocation, Float padding=1.0, Bool usezero=True, 00142 Bool useDoublePrec=False); 00143 GridFT(Long cachesize, Int tilesize, String convType, 00144 MDirection mTangent, Float padding=1.0, Bool usezero=True, 00145 Bool useDoublePrec=False); 00146 GridFT(Long cachesize, Int tilesize, String convType, 00147 MPosition mLocation, MDirection mTangent, Float passing=1.0, 00148 Bool usezero=True, Bool useDoublePrec=False); 00149 // </group> 00150 00151 // Construct from a Record containing the GridFT state 00152 GridFT(const RecordInterface& stateRec); 00153 00154 // Copy constructor 00155 GridFT(const GridFT &other); 00156 00157 // Assignment operator 00158 virtual GridFT &operator=(const GridFT &other); 00159 00160 virtual ~GridFT(); 00161 00162 virtual FTMachine* cloneFTM(); 00163 00164 // Initialize transform to Visibility plane using the image 00165 // as a template. The image is loaded and Fourier transformed. 00166 00167 virtual void initializeToVis(ImageInterface<Complex>& image, 00168 const vi::VisBuffer2& vb); 00169 00170 // Finalize transform to Visibility plane: flushes the image 00171 // cache and shows statistics if it is being used. 00172 virtual void finalizeToVis(); 00173 00174 // Initialize transform to Sky plane: initializes the image 00175 00176 virtual void initializeToSky(ImageInterface<Complex>& image, Matrix<Float>& weight, 00177 const vi::VisBuffer2& vb); 00178 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 virtual void finalizeToSky(); 00183 00184 00185 // Get actual coherence from grid by degridding 00186 00187 virtual void get(vi::VisBuffer2& vb, Int row=-1); 00188 00189 // Put coherence to grid by gridding. 00190 00191 virtual void put(const vi::VisBuffer2& vb, Int row=-1, Bool dopsf=False, 00192 FTMachine::Type type=FTMachine::OBSERVED); 00193 00194 // Make the entire image 00195 void makeImage(FTMachine::Type type, 00196 vi::VisibilityIterator2& vi, 00197 ImageInterface<Complex>& image, 00198 Matrix<Float>& weight); 00199 00200 // Get the final image: do the Fourier transform and 00201 // grid-correct, then optionally normalize by the summed weights 00202 ImageInterface<Complex>& getImage(Matrix<Float>&, Bool normalize=True); 00203 virtual void normalizeImage(Lattice<Complex>& /*skyImage*/, 00204 const Matrix<Double>& /*sumOfWts*/, 00205 Lattice<Float>& /*sensitivityImage*/, 00206 Bool /*fftNorm*/) 00207 {throw(AipsError("GridFT::normalizeImage() called"));} 00208 00209 // Get the final weights image 00210 void getWeightImage(ImageInterface<Float>&, Matrix<Float>&); 00211 00212 // Save and restore the GridFT to and from a record 00213 virtual Bool toRecord(String& error, RecordInterface& outRec, 00214 Bool withImage=False, const String diskimage=""); 00215 virtual Bool fromRecord(String& error, const RecordInterface& inRec); 00216 00217 // Can this FTMachine be represented by Fourier convolutions? 00218 virtual Bool isFourier() {return True;} 00219 00220 virtual void setNoPadding(Bool nopad){noPadding_p=nopad;}; 00221 virtual void modifyConvFunc(const Vector<Double>& convFunc, Int convSupport, Int convSampling); 00222 virtual String name() const; 00223 virtual void setMiscInfo(const Int qualifier){(void)qualifier;}; 00224 virtual void ComputeResiduals(vi::VisBuffer2&/*vb*/, Bool /*useCorrected*/) {}; 00225 00226 protected: 00227 00228 00229 // Padding in FFT 00230 Float padding_p; 00231 00232 // Get the appropriate data pointer 00233 Array<Complex>* getDataPointer(const IPosition&, Bool); 00234 00235 virtual void ok(); 00236 00237 virtual void init(); 00238 00239 //Prepare the grid for degridding 00240 virtual void prepGridForDegrid(); 00241 00242 // Is this record on Grid? check both ends. This assumes that the 00243 // ends bracket the middle 00244 // Bool recordOnGrid(const VisBuffer& vb, Int rownr) const; 00245 00246 00247 // Image cache 00248 LatticeCache<Complex> * imageCache; 00249 00250 // Sizes 00251 Long cachesize; 00252 Int tilesize; 00253 00254 // Gridder 00255 ConvolveGridder<Double, Complex>* gridder; 00256 00257 // Is this tiled? 00258 Bool isTiled; 00259 00260 // Array lattice 00261 SHARED_PTR<Lattice<Complex> > arrayLattice; 00262 00263 // Lattice. For non-tiled gridding, this will point to arrayLattice, 00264 // whereas for tiled gridding, this points to the image 00265 SHARED_PTR<Lattice<Complex> > lattice; 00266 00267 String convType; 00268 00269 Float maxAbsData; 00270 00271 // Useful IPositions 00272 IPosition centerLoc, offsetLoc; 00273 00274 // Image Scaling and offset 00275 Vector<Double> uvScale, uvOffset; 00276 00277 00278 Int priorCacheSize; 00279 00280 // Grid/degrid zero spacing points? 00281 00282 Bool usezero_p; 00283 00284 //force no padding 00285 Bool noPadding_p; 00286 00287 //Check if using put that avoids non-necessary reads 00288 Bool usePut2_p; 00289 00290 //machine name 00291 String machineName_p; 00292 00293 Double timemass_p, timegrid_p, timedegrid_p; 00294 Vector<Double> convFunc_p; 00295 Int convSampling_p, convSupport_p; 00296 // casa::async::SynthesisAsyncPeek *peek; 00297 00298 }; 00299 00300 }//# end of namespace refim 00301 } //# NAMESPACE CASA - END 00302 00303 #endif