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