SDGrid.h

Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1