casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
WOnlyProjectFT.h
Go to the documentation of this file.
00001 //# WOnlyProjectFT.h: Definition for WOnlyProjectFT
00002 //# Copyright (C) 1996,1997,1998,1999,2000,2002
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_WONLYPROJECTFT_H
00030 #define SYNTHESIS_WONLYPROJECTFT_H
00031 
00032 #include <synthesis/MeasurementComponents/VLACalcIlluminationConvFunc.h>
00033 #include <synthesis/MeasurementComponents/VLAIlluminationConvFunc.h>
00034 //#include <synthesis/TransformMachines/ConvolutionFunction.h>
00035 #include <synthesis/MeasurementComponents/EVLAConvFunc.h>
00036 #include <synthesis/MeasurementComponents/SolvableVisCal.h>
00037 #include <synthesis/TransformMachines/VPSkyJones.h>
00038 #include <synthesis/TransformMachines/FTMachine.h>
00039 //#include <synthesis/TransformMachines/CFCache.h>
00040 #include <synthesis/TransformMachines/Utils.h>
00041 
00042 #include <scimath/Mathematics/FFTServer.h>
00043 #include <synthesis/MSVis/VisBuffer.h>
00044 
00045 #include <casa/Containers/Block.h>
00046 #include <casa/Arrays/Array.h>
00047 #include <casa/Arrays/Vector.h>
00048 #include <casa/Arrays/Matrix.h>
00049 
00050 #include <scimath/Mathematics/ConvolveGridder.h>
00051 #include <lattices/Lattices/LatticeCache.h>
00052 #include <lattices/Lattices/ArrayLattice.h>
00053 #include <ms/MeasurementSets/MSColumns.h>
00054 #include <measures/Measures/Measure.h>
00055 #include <measures/Measures/MDirection.h>
00056 #include <measures/Measures/MPosition.h>
00057 #include <images/Images/ImageInterface.h>
00058 #include <coordinates/Coordinates/DirectionCoordinate.h>
00059 
00060 namespace casa { //# NAMESPACE CASA - BEGIN
00061   
00062   // <summary>  An FTMachine for Gridded Fourier transforms including effects of primary beam and pointing offsets and the w-term</summary>
00063   
00064   // <use visibility=export>
00065   
00066   // <reviewed reviewer="" date="" tests="" demos="">
00067   
00068   // <prerequisite>
00069   //   <li> <linkto class=FTMachine>FTMachine</linkto> module
00070   //   <li> <linkto class=SkyEquation>SkyEquation</linkto> module
00071   //   <li> <linkto class=VisBuffer>VisBuffer</linkto> module
00072   //   <li> <linto class=EPJones>EPJones</linkto> module
00073   // </prerequisite>
00074   //
00075   // <etymology> 
00076   // FTMachine is a Machine for Fourier Transforms. Like
00077   // WProjectFT, WOnlyProjectFT does Grid-based Fourier transforms but
00078   // also includes the effects of primary beam and antenna pointing
00079   // offsets. 
00080   // </etymology>
00081   //
00082   // <synopsis> 
00083   //
00084   // The <linkto class=SkyEquation>SkyEquation</linkto> needs to be
00085   // able to perform Fourier transforms on visibility
00086   // data. WOnlyProjectFT allows efficient handling of direction
00087   // dependent effects due to the primary beam and antenna pointing
00088   // offsets using a <linkto class=VisBuffer>VisBuffer</linkto> which
00089   // encapsulates a chunk of visibility (typically all baselines for
00090   // one time) together with all the information needed for processing
00091   // (e.g. UVW coordinates).
00092   //
00093   // Using this FTMachine, errors due antenna pointing offsets can be
00094   // corrected during deconvolution.  One form of antenna pointing
00095   // error which is known a-priori is the VLA polarization squint
00096   // (about 6% of the Primary beam width at any frequency).  For
00097   // Stokes imaging, using this FTMachine, the VLA polarization squint
00098   // and beam polarization can also be corrected.  Also since the
00099   // effects of antenna pointing errors is strongest in the range of
00100   // 1-2GHz band (where the sky is not quite empty while the beams are
00101   // not too large either), this FTMachine can also be setup to
00102   // correct for the w-term.
00103   //
00104   // Switches are provided in the get() method to compute the
00105   // derivatives with respect to the parameters of the primary beam
00106   // (only pointing offsets for now).  This is used in the pointing
00107   // offset solver.
00108   //
00109   // See the documentation of other FTMachines for details about the
00110   // design of the FTMachines in general.
00111   //
00112   // </synopsis> 
00113   //
00114   // <example>
00115   // See the example for <linkto class=SkyModel>SkyModel</linkto>.
00116   // </example>
00117   //
00118   // <motivation>
00119   //
00120   // Encapsulate the correction of direction dependent effects via
00121   // visibility plane convolutions with a potentially different
00122   // convolution function for each baseline.
00123   //
00124   // </motivation>
00125   //
00126   // <todo asof="2005/07/21">
00127   //
00128   // <ul> Include the antenna based complex gain term as well since
00129   // that can interfere with the effects of pointing offsets.  
00130   //
00131   // <ul> Factor out the actual convolution functions as a separate
00132   // class making FTMachines for various direction dependent effects
00133   // generic.
00134   //
00135   // </todo>
00136   
00137   //  class EPJones;
00138   class SolvableVisJones;
00139   class WOnlyProjectFT : public FTMachine {
00140   public:
00141     
00142     // Constructor: cachesize is the size of the cache in words
00143     // (e.g. a few million is a good number), tilesize is the
00144     // size of the tile used in gridding (cannot be less than
00145     // 12, 16 works in most cases). 
00146     // <group>
00147     WOnlyProjectFT(Int nFacets, Long cachesize, 
00148                    CountedPtr<CFCache>& cfcache,
00149                    CountedPtr<ConvolutionFunction>& cf,
00150                    CountedPtr<VisibilityResampler>& reSampler,
00151                    Int tilesize=16, 
00152                    Float pbLimit=5e-2,
00153                    Bool usezero=False);
00154     // </group>
00155     
00156     // Construct from a Record containing the WOnlyProjectFT state
00157     WOnlyProjectFT(const RecordInterface& stateRec);
00158     
00159     // Copy constructor
00160     WOnlyProjectFT(const WOnlyProjectFT &other);
00161     
00162     // Assignment operator
00163     WOnlyProjectFT &operator=(const WOnlyProjectFT &other);
00164     
00165     ~WOnlyProjectFT();
00166     
00167     virtual void initializeToVis(ImageInterface<Complex>& image,
00168                          const VisBuffer& vb);
00169     // This version returns the gridded vis...should be used in conjunction 
00170     // with the version of 'get' that needs the gridded visdata 
00171     virtual void initializeToVis(ImageInterface<Complex>& image,
00172                          const VisBuffer& vb, Array<Complex>& griddedVis,
00173                          Vector<Double>& uvscale);
00174     
00175     // Finalize transform to Visibility plane: flushes the image
00176     // cache and shows statistics if it is being used.
00177     virtual void finalizeToVis();
00178     
00179     // Initialize transform to Sky plane: initializes the image
00180     virtual void initializeToSky(ImageInterface<Complex>& image,  Matrix<Float>& weight,
00181                          const VisBuffer& vb);
00182     
00183     // Finalize transform to Sky plane: flushes the image
00184     // cache and shows statistics if it is being used. DOES NOT
00185     // DO THE FINAL TRANSFORM!
00186     virtual void finalizeToSky();
00187     
00188     virtual void initVisBuffer(VisBuffer& vb, Type whichVBColumn);
00189     void initVisBuffer(VisBuffer& vb, Type whichVBColumn, Int row);
00190 
00191     // Get actual coherence from grid by degridding
00192     void get(VisBuffer& vb, Int row=-1);
00193     
00194     // Get the coherence from grid return it in the degrid 
00195     // is used especially when scratch columns are not 
00196     // present in ms.
00197     void get(VisBuffer& vb, Cube<Complex>& degrid, 
00198              Array<Complex>& griddedVis, Vector<Double>& scale, 
00199              Int row=-1);
00200     
00201     
00202     
00203     // Put coherence to grid by gridding.
00204     void put(const VisBuffer&, TempImage<Complex>&, Vector<Double>&, int,
00205              UVWMachine*, Bool) 
00206     {
00207       //    throw(AipsError("WOnlyProjectFT::put is not implemented"));
00208     }
00209     void put(const VisBuffer& vb, Int row=-1, Bool dopsf=False,
00210              FTMachine::Type type=FTMachine::OBSERVED);
00211     
00212     // Make the entire image
00213     void makeImage(FTMachine::Type type,
00214                    VisSet& vs,
00215                    ImageInterface<Complex>& image,
00216                    Matrix<Float>& weight);
00217     
00218     // Get the final image: do the Fourier transform and
00219     // grid-correct, then optionally normalize by the summed weights
00220     virtual ImageInterface<Complex>& getImage(Matrix<Float>&, Bool normalize=True);
00221 
00222     // Save and restore the WOnlyProjectFT to and from a record
00223     Bool toRecord(RecordInterface& outRec,  Bool withImage=False);
00224     Bool fromRecord(const RecordInterface& inRec);
00225     
00226     // Can this FTMachine be represented by Fourier convolutions?
00227     Bool isFourier() {return True;}
00228     
00229     //
00230     // Make a sensitivity image (sensitivityImage), given the gridded
00231     // weights (wtImage).  These are related to each other by a
00232     // Fourier transform and normalization by the sum-of-weights
00233     // (sumWt) and normalization by the product of the 2D FFT size
00234     // along each axis.  If doFFTNorm=False, normalization by the FFT
00235     // size is not done.  If sumWt is not provided, normalization by
00236     // the sum of weights is also not done.
00237     //
00238     virtual void makeSensitivityImage(Lattice<Complex>& wtImage,
00239                                       ImageInterface<Float>& sensitivityImage,
00240                                       const Matrix<Float>& sumWt=Matrix<Float>(),
00241                                       const Bool& doFFTNorm=True) {};
00242     virtual void makeSensitivityImage(const VisBuffer& vb, 
00243                                       const ImageInterface<Complex>& imageTemplate,
00244                                       ImageInterface<Float>& sensitivityImage);
00245     //
00246     // Given the sky image (Fourier transform of the visibilities),
00247     // sum of weights and the sensitivity image, this method replaces
00248     // the skyImage with the normalized image of the sky.
00249     //
00250     virtual void normalizeImage(Lattice<Complex>& skyImage,
00251                                 const Matrix<Double>& sumOfWts,
00252                                 Lattice<Float>& sensitivityImage,
00253                                 Bool fftNorm=True);
00254     virtual void normalizeImage(Lattice<Complex>& skyImage,
00255                                 const Matrix<Double>& sumOfWts,
00256                                 Lattice<Float>& sensitivityImage,
00257                                 Lattice<Complex>& sensitivitySqImage,
00258                                 Bool fftNorm=True);
00259     
00260     virtual ImageInterface<Float>& getSensitivityImage() {return *avgPB_p;}
00261     virtual Matrix<Double>& getSumOfWeights() {return sumWeight;};
00262 
00263     Vector<Int>& getPolMap() {return polMap;};
00264     virtual String name(){ return "WOnlyProjectFT";};
00265     virtual Bool verifyShapes(IPosition shape0, IPosition shape1);
00266     virtual void setMiscInfo(const Int qualifier){(void)qualifier;};
00267 
00268   protected:
00269     
00270     Int nint(Double val) {return Int(floor(val+0.5));};
00271     // Locate convolution functions on the disk
00272     //    Int locateConvFunction(const Int Nw, const Float pa);
00273     //    void cacheConvFunction(Int which, Array<Complex>& cf, CoordinateSystem& coord);
00274     // Find the convolution function
00275     void findConvFunction(const ImageInterface<Complex>& image,
00276                           const VisBuffer& vb);
00277     
00278     // Get the appropriate data pointer
00279     Array<Complex>* getDataPointer(const IPosition&, Bool);
00280     
00281     void ok();
00282     
00283     void init();
00284     //    virtual void initPolInfo(const VisBuffer& vb);
00285     // Is this record on Grid? check both ends. This assumes that the
00286     // ends bracket the middle
00287     Bool recordOnGrid(const VisBuffer& vb, Int rownr) const;
00288     
00289     // Padding in FFT
00290     Float padding_p;
00291     
00292     Int nWPlanes_p;
00293     // Image cache
00294     LatticeCache<Complex> * imageCache;
00295     
00296     // Sizes
00297     Long cachesize;
00298     Int tilesize;
00299     
00300     // Gridder
00301     ConvolveGridder<Double, Complex>* gridder;
00302     
00303     // Is this tiled?
00304     Bool isTiled;
00305     
00306     // Array lattice
00307     CountedPtr<Lattice<Complex> > arrayLattice;
00308     
00309     // Lattice. For non-tiled gridding, this will point to arrayLattice,
00310     //  whereas for tiled gridding, this points to the image
00311     CountedPtr<Lattice<Complex> > lattice;
00312     
00313     Float maxAbsData;
00314     
00315     // Useful IPositions
00316     IPosition centerLoc, offsetLoc;
00317     
00318     // Image Scaling and offset
00319     Vector<Double> uvScale, uvOffset;
00320     
00321     // Array for non-tiled gridding
00322     Array<Complex> griddedData;
00323     
00324     //    DirectionCoordinate directionCoord;
00325     MDirection::Convert* pointingToImage;
00326     
00327     // Grid/degrid zero spacing points?
00328     Bool usezero_p;
00329     
00330     //    CountedPtr<ConvolutionFunction> telescopeConvFunc_p;
00331     //    CFStore cfs_p, cfwts_p;
00332     Array<Complex> convFunc_p, convWeights_p;
00333     CountedPtr<VisibilityResampler> reSampler_p;
00334     //
00335     // Vector to hold the support size info. for the convolution
00336     // functions pointed to by the elements of convFunctions_p.  The
00337     // co-ordinates of this array are (W-term, Poln, PA).
00338     //
00339     Int convSize, convSampling, wConvSize, lastIndex_p;
00340     
00341     //
00342     // No. of vis. polarization planes used in making the user defined
00343     // Stokes images
00344     //
00345     Int maxConvSupport;
00346 
00347     Int Nant_p;
00348     Bool makingPSF;
00349     
00350     virtual void runFortranGet(Matrix<Double>& uvw,Vector<Double>& dphase,
00351                                Cube<Complex>& visdata,
00352                                IPosition& s,
00353                                //Cube<Complex>& gradVisAzData,
00354                                //Cube<Complex>& gradVisElData,
00355                                //IPosition& gradS,
00356                                Int& Conj,
00357                                Cube<Int>& flags,Vector<Int>& rowFlags,
00358                                Int& rownr,Vector<Double>& actualOffset,
00359                                Array<Complex>* dataPtr,
00360                                Int& aNx, Int& aNy, Int& npol, Int& nchan,
00361                                VisBuffer& vb,Int& Nant_p, Int& scanNo,
00362                                Double& sigma,
00363                                Array<Float>& raoffsets,
00364                                Array<Float>& decoffsets,
00365                                Double area,
00366                                Int& doGrad,Int paIndex);
00367     virtual void runFortranPut(Matrix<Double>& uvw,Vector<Double>& dphase,
00368                                const Complex& visdata_p,
00369                                IPosition& s,
00370                                //Cube<Complex>& gradVisAzData,
00371                                //Cube<Complex>& gradVisElData,
00372                                //IPosition& gradS,
00373                        Int& Conj,
00374                        Cube<Int>& flags,Vector<Int>& rowFlags,
00375                        const Matrix<Float>& weight,
00376                        Int& rownr,Vector<Double>& actualOffset,
00377                        Array<Complex>& dataPtr,
00378                        Int& aNx, Int& aNy, Int& npol, Int& nchan,
00379                        const VisBuffer& vb,Int& Nant_p, Int& scanNo,
00380                        Double& sigma,
00381                        Array<Float>& raoffsets,
00382                        Array<Float>& decoffsets,
00383                        Matrix<Double>& sumWeight,
00384                        Double& area,
00385                        Int& doGrad,
00386                        Int& doPSF,Int paIndex);
00387   void runFortranGetGrad(Matrix<Double>& uvw,Vector<Double>& dphase,
00388                          Cube<Complex>& visdata,
00389                          IPosition& s,
00390                          Cube<Complex>& gradVisAzData,
00391                          Cube<Complex>& gradVisElData,
00392                          //                      IPosition& gradS,
00393                          Int& Conj,
00394                          Cube<Int>& flags,Vector<Int>& rowFlags,
00395                          Int& rownr,Vector<Double>& actualOffset,
00396                          Array<Complex>* dataPtr,
00397                          Int& aNx, Int& aNy, Int& npol, Int& nchan,
00398                          VisBuffer& vb,Int& Nant_p, Int& scanNo,
00399                          Double& sigma,
00400                          Array<Float>& l_off,
00401                          Array<Float>& m_off,
00402                          Double area,
00403                          Int& doGrad,
00404                          Int paIndex);
00405   };
00406 } //# NAMESPACE CASA - END
00407     // void makeAntiAliasingOp(Vector<Complex>& val, const Int len, const Double HPBW);
00408     // void makeAntiAliasingCorrection(Vector<Complex>& correction, 
00409     //                              const Vector<Complex>& op, 
00410     //                              const Int nx);
00411     // void applyAntiAliasingOp(ImageInterface<Complex>& cf, 
00412     //                       Vector<IPosition>& offset,
00413     //                       Double HPBW,
00414     //                       Int op=0, 
00415     //                       Bool Square=False);
00416     // void applyAntiAliasingOp(ImageInterface<Float>& cf, 
00417     //                       Vector<IPosition>& offset,
00418     //                       Double HPBW,
00419     //                       Int op=0, 
00420     //                       Bool Square=False);
00421     // void correctAntiAliasing(Lattice<Complex>& cf);
00422 
00423 #endif