00001 //# SDGrid.h: Definition for SDGrid 00002 //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 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_SDGRID_H 00030 #define SYNTHESIS_SDGRID_H 00031 00032 #include <casa/Arrays/Array.h> 00033 #include <casa/Arrays/Matrix.h> 00034 #include <casa/Arrays/Vector.h> 00035 #include <casa/Containers/Block.h> 00036 #include <coordinates/Coordinates/DirectionCoordinate.h> 00037 #include <images/Images/ImageInterface.h> 00038 #include <lattices/Lattices/ArrayLattice.h> 00039 #include <lattices/Lattices/LatticeCache.h> 00040 #include <measures/Measures/Measure.h> 00041 #include <measures/Measures/MDirection.h> 00042 #include <measures/Measures/MPosition.h> 00043 #include <ms/MeasurementSets/MSColumns.h> 00044 #include <msvis/MSVis/VisBuffer.h> 00045 #include <scimath/Mathematics/FFTServer.h> 00046 #include <synthesis/MeasurementComponents/SDPosInterpolator.h> 00047 #include <synthesis/TransformMachines/FTMachine.h> 00048 #include <synthesis/TransformMachines/SkyJones.h> 00049 00050 namespace casa { //# NAMESPACE CASA - BEGIN 00051 00052 // <summary> An FTMachine for Gridding Single Dish data 00053 // </summary> 00054 00055 // <use visibility=export> 00056 00057 // <reviewed reviewer="" date="" tests="" demos=""> 00058 00059 // <prerequisite> 00060 // <li> <linkto class=FTMachine>FTMachine</linkto> module 00061 // <li> <linkto class=SkyEquation>SkyEquation</linkto> module 00062 // <li> <linkto class=VisBuffer>VisBuffer</linkto> module 00063 // </prerequisite> 00064 // 00065 // <etymology> 00066 // FTMachine is a Machine for Fourier Transforms. SDGrid does 00067 // Single Dish gridding in a similar way 00068 // </etymology> 00069 // 00070 // <synopsis> 00071 // The <linkto class=SkyEquation>SkyEquation</linkto> needs to be able 00072 // to perform Fourier transforms on visibility data and to grid 00073 // single dish data. 00074 // SDGrid allows efficient Single Dish processing using a 00075 // <linkto class=VisBuffer>VisBuffer</linkto> which encapsulates 00076 // a chunk of visibility (typically all baselines for one time) 00077 // together with all the information needed for processing 00078 // (e.g. direction coordinates). 00079 // 00080 // Gridding and degridding in SDGrid are performed using a 00081 // novel sort-less algorithm. In this approach, the gridded plane is 00082 // divided into small patches, a cache of which is maintained in memory 00083 // using a general-purpose <linkto class=LatticeCache>LatticeCache</linkto> class. As the (time-sorted) 00084 // visibility data move around slowly in the image plane, patches are 00085 // swapped in and out as necessary. Thus, optimally, one would keep at 00086 // least one patch per scan line of data. 00087 // 00088 // A grid cache is defined on construction. If the gridded image plane is smaller 00089 // than this, it is kept entirely in memory and all gridding and 00090 // degridding is done entirely in memory. Otherwise a cache of tiles is 00091 // kept an paged in and out as necessary. Optimally the cache should be 00092 // big enough to hold all polarizations and frequencies for one 00093 // complete scan line. 00094 // The paging rate will then be small. As the cache size is 00095 // reduced below this critical value, paging increases. The algorithm will 00096 // work for only one patch but it will be very slow! 00097 // 00098 // The gridding and degridding steps are implemented in Fortran 00099 // for speed. In gridding, the visibilities are added onto the 00100 // grid points in the neighborhood using a weighting function. 00101 // In degridding, the value is derived by a weight summ of the 00102 // same points, using the same weighting function. 00103 // </synopsis> 00104 // 00105 // <example> 00106 // See the example for <linkto class=SkyModel>SkyModel</linkto>. 00107 // </example> 00108 // 00109 // <motivation> 00110 // Define an interface to allow efficient processing of chunks of 00111 // visibility data 00112 // </motivation> 00113 // 00114 // <todo asof="97/10/01"> 00115 // <ul> Deal with large VLA spectral line case 00116 // </todo> 00117 00118 class SDGrid : public FTMachine { 00119 public: 00120 00121 // Constructor: cachesize is the size of the cache in words 00122 // (e.g. a few million is a good number), tilesize is the 00123 // size of the tile used in gridding (cannot be less than 00124 // 12, 16 works in most cases), and convType is the type of 00125 // gridding used (SF is prolate spheriodal wavefunction, 00126 // and BOX is plain box-car summation). mLocation is 00127 // the position to be used in some phase rotations. If 00128 // mTangent is specified then the uvw rotation is done for 00129 // that location iso the image center. userSupport is to allow 00130 // larger support for the convolution if the user wants it ..-1 will 00131 // use the default i.e 1 for BOX and 3 for others 00132 // USEIMAGINGWEIGHT 00133 // The parameter useImagingWeight in the constructors is to explicitly 00134 // use vb.imagingweight while gridding, 00135 // When doing just SD imaging then setting it to False is fine (in fact recommended as vb.imagingweight 00136 // is set to zero if any pol is flagged this may change later .....today being 2014/08/06) 00137 // when using it in conjuction with interferometer gridding then set useImagingWeight to True 00138 // this is to allow for proper non natural weighting scheme while imaging 00139 // <group> 00140 SDGrid(SkyJones& sj, Int cachesize, Int tilesize, 00141 String convType="BOX", Int userSupport=-1, Bool useImagingWeight=False); 00142 SDGrid(MPosition& ml, SkyJones& sj, Int cachesize, 00143 Int tilesize, String convType="BOX", Int userSupport=-1, 00144 Float minweight=0., Bool clipminmax=False, Bool useImagingWeight=False); 00145 SDGrid(Int cachesize, Int tilesize, 00146 String convType="BOX", Int userSupport=-1, Bool useImagingWeight=False); 00147 SDGrid(MPosition& ml, Int cachesize, Int tilesize, 00148 String convType="BOX", Int userSupport=-1, Float minweight=0., Bool clipminmax=False, 00149 Bool useImagingWeight=False); 00150 SDGrid(MPosition& ml, Int cachesize, Int tilesize, 00151 String convType="TGAUSS", Float truncate=-1.0, 00152 Float gwidth=0.0, Float jwidth=0.0, Float minweight=0., Bool clipminmax=False, 00153 Bool useImagingWeight=False); 00154 // </group> 00155 00156 // Copy constructor 00157 SDGrid(const SDGrid &other); 00158 00159 // Assignment operator 00160 SDGrid &operator=(const SDGrid &other); 00161 00162 ~SDGrid(); 00163 00164 // Initialize transform to Visibility plane using the image 00165 // as a template. The image is loaded and Fourier transformed. 00166 void initializeToVis(ImageInterface<Complex>& image, 00167 const VisBuffer& vb); 00168 00169 // Finalize transform to Visibility plane: flushes the image 00170 // cache and shows statistics if it is being used. 00171 void finalizeToVis(); 00172 00173 // Initialize transform to Sky plane: initializes the image 00174 void initializeToSky(ImageInterface<Complex>& image, Matrix<Float>& weight, 00175 const VisBuffer& vb); 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 void finalizeToSky(); 00181 00182 // Get actual coherence from grid by degridding 00183 void get(VisBuffer& vb, Int row=-1); 00184 00185 // Put coherence to grid by gridding. 00186 void put(const VisBuffer& vb, Int row=-1, Bool dopsf=False, 00187 FTMachine::Type type=FTMachine::OBSERVED); 00188 00189 // Make the entire image using a ROVisIter... 00190 // This is an overload for FTMachine version as 00191 //SDGrid now does everything in memory 00192 // so for large cube ..proceed by slices that fit in memory here. 00193 virtual void makeImage(FTMachine::Type type, 00194 ROVisibilityIterator& vi, 00195 ImageInterface<Complex>& image, 00196 Matrix<Float>& weight); 00197 00198 // Get the final image: 00199 // optionally normalize by the summed weights 00200 ImageInterface<Complex>& getImage(Matrix<Float>&, Bool normalize=True); 00201 virtual void normalizeImage(Lattice<Complex>& /*skyImage*/, 00202 const Matrix<Double>& /*sumOfWts*/, 00203 Lattice<Float>& /*sensitivityImage*/, 00204 Bool /*fftNorm*/) 00205 {throw(AipsError("SDGrid::normalizeImage() called"));} 00206 00207 // Get the final weights image 00208 void getWeightImage(ImageInterface<Float>&, Matrix<Float>&); 00209 00210 // Has this operator changed since the last application? 00211 virtual Bool changed(const VisBuffer& vb); 00212 virtual void setMiscInfo(const Int qualifier){(void)qualifier;}; 00213 virtual void ComputeResiduals(VisBuffer& /*vb*/, Bool /*useCorrected*/) {}; 00214 00215 virtual String name() const; 00216 00217 private: 00218 00219 // Find the Primary beam and convert it into a convolution buffer 00220 void findPBAsConvFunction(const ImageInterface<Complex>& image, 00221 const VisBuffer& vb); 00222 00223 SkyJones* sj_p; 00224 00225 // Get the appropriate data pointer 00226 Array<Complex>* getDataPointer(const IPosition&, Bool); 00227 Array<Float>* getWDataPointer(const IPosition&, Bool); 00228 00229 void ok(); 00230 00231 void init(); 00232 00233 // Image cache 00234 LatticeCache<Complex> * imageCache; 00235 LatticeCache<Float> * wImageCache; 00236 00237 // Sizes 00238 Int cachesize, tilesize; 00239 00240 // Is this tiled? 00241 Bool isTiled; 00242 00243 // Storage for weights 00244 ImageInterface<Float>* wImage; 00245 00246 // Array lattice 00247 Lattice<Complex> * arrayLattice; 00248 Lattice<Float> * wArrayLattice; 00249 00250 // Lattice. For non-tiled gridding, this will point to arrayLattice, 00251 // whereas for tiled gridding, this points to the image 00252 Lattice<Complex>* lattice; 00253 Lattice<Float>* wLattice; 00254 00255 String convType; 00256 00257 // Useful IPositions 00258 IPosition centerLoc, offsetLoc; 00259 00260 // Array for non-tiled gridding 00261 Array<Float> wGriddedData; 00262 00263 00264 DirectionCoordinate directionCoord; 00265 00266 MDirection::Convert* pointingToImage; 00267 00268 Vector<Double> xyPos; 00269 //Original xypos of moving source 00270 Vector<Double> xyPosMovingOrig_p; 00271 00272 MDirection worldPosMeas; 00273 00274 Cube<Int> flags; 00275 00276 Vector<Float> convFunc; 00277 Int convSampling; 00278 Int convSize; 00279 Int convSupport; 00280 Int userSetSupport_p; 00281 00282 Float truncate_p; 00283 Float gwidth_p; 00284 Float jwidth_p; 00285 00286 Float minWeight_p; 00287 00288 Int lastIndex_p; 00289 Bool useImagingWeight_p; 00290 Int lastAntID_p; 00291 Int msId_p; 00292 00293 Bool isSplineInterpolationReady; 00294 SDPosInterpolator* interpolator; 00295 00296 // for minmax clipping 00297 Bool clipminmax_; 00298 Array<Complex> gmin_; 00299 Array<Complex> gmax_; 00300 Array<Float> wmin_; 00301 Array<Float> wmax_; 00302 Array<Int> npoints_; 00303 void clipMinMax(); 00304 00305 Int getIndex(const ROMSPointingColumns& mspc, const Double& time, 00306 const Double& interval=-1.0, const Int& antid=-1); 00307 00308 Bool getXYPos(const VisBuffer& vb, Int row); 00309 00310 //get the MDirection from a chosen column of pointing table 00311 MDirection directionMeas(const ROMSPointingColumns& mspc, const Int& index); 00312 MDirection directionMeas(const ROMSPointingColumns& mspc, const Int& index, const Double& time); 00313 MDirection interpolateDirectionMeas(const ROMSPointingColumns& mspc, const Double& time, 00314 const Int& index, const Int& index1, const Int& index2); 00315 00316 void pickWeights(const VisBuffer&vb, Matrix<Float>& weight); 00317 00318 //for debugging 00319 //FILE *pfile; 00320 }; 00321 00322 } //# NAMESPACE CASA - END 00323 00324 #endif