casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
AWProjectFT.h
Go to the documentation of this file.
00001 //# AWProjectFT.h: Definition for AWProjectFT
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_AWPROJECTFT_H
00030 #define SYNTHESIS_AWPROJECTFT_H
00031 
00032 #include <synthesis/TransformMachines/VLACalcIlluminationConvFunc.h>
00033 #include <synthesis/TransformMachines/VLAIlluminationConvFunc.h>
00034 #include <synthesis/TransformMachines/AWVisResampler.h>
00035 //#include <synthesis/MeasurementComponents/ConvolutionFunction.h>
00036 #include <synthesis/TransformMachines/EVLAConvFunc.h>
00037 #include <synthesis/MeasurementComponents/SolvableVisCal.h>
00038 #include <synthesis/TransformMachines/VPSkyJones.h>
00039 #include <synthesis/TransformMachines/FTMachine.h>
00040 //#include <synthesis/MeasurementComponents/CFCache.h>
00041 #include <synthesis/TransformMachines/Utils.h>
00042 
00043 #include <scimath/Mathematics/FFTServer.h>
00044 #include <synthesis/MSVis/VisBuffer.h>
00045 
00046 #include <casa/Containers/Block.h>
00047 #include <casa/Arrays/Array.h>
00048 #include <casa/Arrays/Vector.h>
00049 #include <casa/Arrays/Matrix.h>
00050 
00051 #include <scimath/Mathematics/ConvolveGridder.h>
00052 #include <lattices/Lattices/LatticeCache.h>
00053 #include <lattices/Lattices/ArrayLattice.h>
00054 #include <ms/MeasurementSets/MSColumns.h>
00055 #include <measures/Measures/Measure.h>
00056 #include <measures/Measures/MDirection.h>
00057 #include <measures/Measures/MPosition.h>
00058 #include <images/Images/ImageInterface.h>
00059 #include <coordinates/Coordinates/DirectionCoordinate.h>
00060 
00061 #include <casa/OS/Timer.h>
00062 
00063 namespace casa { //# NAMESPACE CASA - BEGIN
00064   
00065   // <summary>  An FTMachine for Gridded Fourier transforms including effects of primary beam and pointing offsets and the w-term</summary>
00066   
00067   // <use visibility=export>
00068   
00069   // <reviewed reviewer="" date="" tests="" demos="">
00070   
00071   // <prerequisite>
00072   //   <li> <linkto class=FTMachine>FTMachine</linkto> module
00073   //   <li> <linkto class=SkyEquation>SkyEquation</linkto> module
00074   //   <li> <linkto class=VisBuffer>VisBuffer</linkto> module
00075   //   <li> <linto class=EPJones>EPJones</linkto> module
00076   // </prerequisite>
00077   //
00078   // <etymology> 
00079   // FTMachine is a Machine for Fourier Transforms. Like
00080   // WProjectFT, AWProjectFT does Grid-based Fourier transforms but
00081   // also includes the effects of primary beam and antenna pointing
00082   // offsets. 
00083   // </etymology>
00084   //
00085   // <synopsis> 
00086   //
00087   // The <linkto class=SkyEquation>SkyEquation</linkto> needs to be
00088   // able to perform Fourier transforms on visibility
00089   // data. AWProjectFT allows efficient handling of direction
00090   // dependent effects due to the primary beam and antenna pointing
00091   // offsets using a <linkto class=VisBuffer>VisBuffer</linkto> which
00092   // encapsulates a chunk of visibility (typically all baselines for
00093   // one time) together with all the information needed for processing
00094   // (e.g. UVW coordinates).
00095   //
00096   // Using this FTMachine, errors due antenna pointing offsets can be
00097   // corrected during deconvolution.  One form of antenna pointing
00098   // error which is known a-priori is the VLA polarization squint
00099   // (about 6% of the Primary beam width at any frequency).  For
00100   // Stokes imaging, using this FTMachine, the VLA polarization squint
00101   // and beam polarization can also be corrected.  Also since the
00102   // effects of antenna pointing errors is strongest in the range of
00103   // 1-2GHz band (where the sky is not quite empty while the beams are
00104   // not too large either), this FTMachine can also be setup to
00105   // correct for the w-term.
00106   //
00107   // Switches are provided in the get() method to compute the
00108   // derivatives with respect to the parameters of the primary beam
00109   // (only pointing offsets for now).  This is used in the pointing
00110   // offset solver.
00111   //
00112   // See the documentation of other FTMachines for details about the
00113   // design of the FTMachines in general.
00114   //
00115   // </synopsis> 
00116   //
00117   // <example>
00118   // See the example for <linkto class=SkyModel>SkyModel</linkto>.
00119   // </example>
00120   //
00121   // <motivation>
00122   //
00123   // Encapsulate the correction of direction dependent effects via
00124   // visibility plane convolutions with a potentially different
00125   // convolution function for each baseline.
00126   //
00127   // </motivation>
00128   //
00129   // <todo asof="2005/07/21">
00130   //
00131   // <ul> Include the antenna based complex gain term as well since
00132   // that can interfere with the effects of pointing offsets.  
00133   //
00134   // <ul> Factor out the actual convolution functions as a separate
00135   // class making FTMachines for various direction dependent effects
00136   // generic.
00137   //
00138   // </todo>
00139   
00140   //  class EPJones;
00141   class SolvableVisJones;
00142   class AWProjectFT : public FTMachine {
00143   public:
00144     AWProjectFT();
00145     
00146     // Constructor: cachesize is the size of the cache in words
00147     // (e.g. a few million is a good number), tilesize is the
00148     // size of the tile used in gridding (cannot be less than
00149     // 12, 16 works in most cases). 
00150     // <group>
00151     AWProjectFT(Int nFacets, Long cachesize, 
00152                 CountedPtr<CFCache>& cfcache,
00153                 CountedPtr<ConvolutionFunction>& cf,
00154                 CountedPtr<VisibilityResamplerBase>& visResampler,
00155                 Bool applyPointingOffset=True,
00156                 Bool doPBCorr=True,
00157                 Int tilesize=16, 
00158                 Float pbLimit=5e-4,
00159                 Bool usezero=False,
00160                 Bool conjBeams_p=True,
00161                 Bool doublePrecGrid=False);
00162     // </group>
00163     
00164     // Construct from a Record containing the AWProjectFT state
00165     AWProjectFT(const RecordInterface& stateRec);
00166     
00167     // Copy constructor
00168     AWProjectFT(const AWProjectFT &other);
00169     
00170     // Assignment operator
00171     AWProjectFT &operator=(const AWProjectFT &other);
00172     
00173     ~AWProjectFT();
00174     
00175     //   void setEPJones(EPJones* ep_j) {epJ = ep_j;}
00176     void setEPJones(SolvableVisJones* ep_j) {epJ_p = ep_j;}
00177     
00178     virtual void setDOPBCorrection(Bool doit=True) {doPBCorrection=doit;};
00179     virtual Bool getDOPBCorrection() {return doPBCorrection;};
00180     virtual void setConjBeams(Bool useit=True) {conjBeams_p=useit;};
00181     virtual Bool getConjBeams() {return conjBeams_p;};
00182 
00183     virtual Float getPBLimit() {return pbLimit_p;};
00184     // Initialize transform to Visibility plane using the image
00185     // as a template. The image is loaded and Fourier transformed.
00186     
00187     void setObservatoryLocation(const MPosition& mLocation) {mLocation_p=mLocation;};
00188     
00189     // Vectorized version of initializeToVis.  Required since
00190     // MultiTermFTM needs vectorized version. And this is implemented
00191     // in FTMachine (the baseclass).  And one relies on static_cast in
00192     // Imager::createFTMachine() (or is it in CSE?) to cast the
00193     // pointer to specific types such that this methods gets called
00194     // when the FTMachine pointer is of type AWProjectFT.
00195     virtual void initializeToVis(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec,
00196                                  PtrBlock<SubImage<Float> *> & modelImageVec, 
00197                                  PtrBlock<SubImage<Float> *>& weightImageVec, 
00198                                  PtrBlock<SubImage<Float> *>& fluxScaleVec, 
00199                                  Block<Matrix<Float> >& weightsVec,
00200                                  const VisBuffer& vb);
00201 
00202     virtual void initializeToVis(ImageInterface<Complex>& image,
00203                          const VisBuffer& vb);
00204     // This version returns the gridded vis...should be used in conjunction 
00205     // with the version of 'get' that needs the gridded visdata 
00206     virtual void initializeToVis(ImageInterface<Complex>& image,
00207                          const VisBuffer& vb, Array<Complex>& griddedVis,
00208                          Vector<Double>& uvscale);
00209     
00210     // Finalize transform to Visibility plane: flushes the image
00211     // cache and shows statistics if it is being used.
00212     virtual void finalizeToVis();
00213     
00214     // Initialize transform to Sky plane: initializes the image
00215     virtual void initializeToSky(ImageInterface<Complex>& image,  Matrix<Float>& weight,
00216                          const VisBuffer& vb);
00217     
00218     // Finalize transform to Sky plane: flushes the image
00219     // cache and shows statistics if it is being used. DOES NOT
00220     // DO THE FINAL TRANSFORM!
00221     virtual void finalizeToSky();
00222     
00223     virtual void initVisBuffer(VisBuffer& vb, Type whichVBColumn);
00224     void initVisBuffer(VisBuffer& vb, Type whichVBColumn, Int row);
00225 
00226     // Get actual coherence from grid by degridding
00227     void get(VisBuffer& vb, Int row=-1);
00228     
00229     // Get the coherence from grid return it in the degrid 
00230     // is used especially when scratch columns are not 
00231     // present in ms.
00232     void get(VisBuffer& vb, Cube<Complex>& degrid, 
00233              Array<Complex>& griddedVis, Vector<Double>& scale, 
00234              Int row=-1);
00235     
00236     void get(VisBuffer& vb, Cube<Float>& pointingOffsets, Int row=-1,
00237              Type whichVBColumn=FTMachine::MODEL,Int Conj=0)
00238     {
00239       get(vb,vb,vb,pointingOffsets,row,whichVBColumn,whichVBColumn,Conj,0);
00240     }
00241     
00242     void get(VisBuffer& vb, VisBuffer& gradAzVB,VisBuffer& gradElVB,
00243              Cube<Float>& pointingOffsets,Int row=-1,
00244              Type whichVBColumn=FTMachine::MODEL,
00245              Type whichGradVBColumn=FTMachine::MODEL,
00246              Int Conj=0, Int doGrad=1) ;
00247   void nget(VisBuffer& vb,
00248             // These offsets should be appropriate for the VB
00249             Array<Float>& l_off, Array<Float>& m_off,
00250             Cube<Complex>& Mout,
00251             Cube<Complex>& dMout1,
00252             Cube<Complex>& dMout2,
00253             Int Conj=0, Int doGrad=1);
00254     // Get the coherence from grid return it in the degrid 
00255     // is used especially when scratch columns are not 
00256     // present in ms.
00257     void get(VisBuffer& vb, Cube<Complex>& degrid, 
00258              Array<Complex>& griddedVis, Vector<Double>& scale, 
00259              Cube<Float>& pointingOffsets,Int row=-1);
00260     
00261     
00262     // Put coherence to grid by gridding.
00263     void put(const VisBuffer&,
00264              TempImage<Complex>&, Vector<Double>&, int,
00265              UVWMachine*, Bool) 
00266     {
00267       //    throw(AipsError("AWProjectFT::put is not implemented"));
00268     }
00269     void put(const VisBuffer& vb, Int row=-1, Bool dopsf=False,
00270              FTMachine::Type type=FTMachine::OBSERVED);
00271     
00272     // Make the entire image using a ROVisIter
00273     virtual void makeImage(FTMachine::Type,
00274                            ROVisibilityIterator&,
00275                            ImageInterface<Complex>&,
00276                            Matrix<Float>&) {};
00277 
00278     // Make the entire image
00279     void makeImage(FTMachine::Type type,
00280                    VisSet& vs,
00281                    ImageInterface<Complex>& image,
00282                    Matrix<Float>& weight);
00283     
00284     // Get the final image: do the Fourier transform and
00285     // grid-correct, then optionally normalize by the summed weights
00286     virtual ImageInterface<Complex>& getImage(Matrix<Float>&, Bool normalize=True);
00287 
00288     // Get the final weights image
00289     void getWeightImage(ImageInterface<Float>&, Matrix<Float>&);
00290     
00291     // Save and restore the AWProjectFT to and from a record
00292     Bool toRecord(RecordInterface& outRec,  Bool withImage=False);
00293     Bool fromRecord(const RecordInterface& inRec);
00294     
00295     // Can this FTMachine be represented by Fourier convolutions?
00296     Bool isFourier() {return True;}
00297     
00298     //  Bool changed(const VisBuffer& vb) {return vpSJ->changed(vb,1);};
00299     //    Bool changed(const VisBuffer& vb) {return False;}
00300     
00301     virtual Int findPointingOffsets(const VisBuffer&, Array<Float>&, Array<Float>&,
00302                                     Bool Evaluate=True);
00303     virtual Int findPointingOffsets(const VisBuffer&, Cube<Float>&,
00304                                     Array<Float>&, Array<Float>&,
00305                                     Bool Evaluate=True);
00306     virtual Double getVBPA(const VisBuffer& vb) 
00307     {
00308       // if (!rotateApertureOTP_p) return currentCFPA;
00309       // else return getPA(vb);
00310       return getPA(vb);
00311     };
00312     MDirection::Convert makeCoordinateMachine(const VisBuffer&,
00313                                               const MDirection::Types&,
00314                                               const MDirection::Types&,
00315                                               MEpoch& last);
00316     //
00317     // Make a sensitivity image (sensitivityImage), given the gridded
00318     // weights (wtImage).  These are related to each other by a
00319     // Fourier transform and normalization by the sum-of-weights
00320     // (sumWt) and normalization by the product of the 2D FFT size
00321     // along each axis.  If doFFTNorm=False, normalization by the FFT
00322     // size is not done.  If sumWt is not provided, normalization by
00323     // the sum of weights is also not done.
00324     //
00325     virtual void makeSensitivityImage(Lattice<Complex>&,// wtImage,
00326                                       ImageInterface<Float>&,// sensitivityImage,
00327                                       const Matrix<Float>&,// sumWt=Matrix<Float>(),
00328                                       const Bool& ) {};
00329     virtual void makeSensitivityImage(const VisBuffer& vb, const ImageInterface<Complex>& imageTemplate,
00330                                       ImageInterface<Float>& sensitivityImage);
00331 
00332     //
00333     // Given the sky image (Fourier transform of the visibilities),
00334     // sum of weights and the sensitivity image, this method replaces
00335     // the skyImage with the normalized image of the sky.
00336     //
00337     // These are now implement in the FTMachine base class.  
00338     //
00339     // These can't be left here since now the design depends on
00340     // getting the correct version called due to static type casting
00341     // of FTM pointer and not on inheretance and and specialization
00342     // via overloading.
00343 
00344     // virtual void normalizeImage(Lattice<Complex>& skyImage,
00345     //                          const Matrix<Double>& sumOfWts,
00346     //                          Lattice<Float>& sensitivityImage,
00347     //                          Bool fftNorm=True);
00348     // virtual void normalizeImage(Lattice<Complex>& skyImage,
00349     //                          const Matrix<Double>& sumOfWts,
00350     //                          Lattice<Float>& sensitivityImage,
00351     //                          Lattice<Complex>& sensitivitySqImage,
00352     //                          Bool fftNorm=True);
00353     
00354     virtual ImageInterface<Float>& getSensitivityImage() {return *avgPB_p;}
00355     virtual Matrix<Double>& getSumOfWeights() {return sumWeight;};
00356     virtual Matrix<Double>& getSumOfCFWeights() {return sumCFWeight;};
00357 
00358     void makeConjPolMap(const VisBuffer& vb, const Vector<Int> cfPolMap, Vector<Int>& conjPolMap);
00359     //    Vector<Int> makeConjPolMap(const VisBuffer& vb);
00360     void makeCFPolMap(const VisBuffer& vb, const Vector<Int>& cfstokes, Vector<Int>& polM);
00361     //    void reset() {vpSJ->reset();}
00362     void reset() {paChangeDetector.reset();}
00363 
00364     void setPAIncrement(const Quantity &computePAIncr, const Quantity &rotateOTFPAIncr);
00365 
00366     Vector<Int>& getPolMap() {return polMap;};
00367     virtual String name() const { return "AWProjectFT";};
00368     virtual Bool verifyAvgPB(ImageInterface<Float>& pb, ImageInterface<Float>& sky)
00369     {return verifyShapes(pb.shape(),sky.shape());}
00370 
00371     virtual Bool verifyAvgPB(ImageInterface<Float>& pb, ImageInterface<Complex>& sky)
00372     {return verifyShapes(pb.shape(),sky.shape());}
00373 
00374     virtual Bool verifyShapes(IPosition shape0, IPosition shape1);
00375 
00376     inline virtual Float pbFunc(const Float& a, const Float& limit) 
00377     //    {Float tt=sqrt(a);return  (abs(tt) >= limit)?tt:1.0;};
00378     {if (abs(a) >= limit) return (a);else return 1.0;};
00379     inline virtual Complex pbFunc(const Complex& a, const Float& limit) 
00380     {if (abs(a)>=limit) return (a); else return Complex(1.0,0.0);};
00381 
00382     virtual void setMiscInfo(const Int qualifier)
00383     {
00384       sensitivityPatternQualifier_p=qualifier;
00385       sensitivityPatternQualifierStr_p = "_tt"+String::toString(sensitivityPatternQualifier_p);
00386     }
00387     virtual void ComputeResiduals(VisBuffer&vb, Bool useCorrected);
00388     void makeWBCFWt(CFStore2& cfs,const Double imRefFreq);
00389 
00390   protected:
00391     
00392     Int nint(Double val) {return Int(floor(val+0.5));};
00393     // Locate convolution functions on the disk
00394     //    Int locateConvFunction(const Int Nw, const Float pa);
00395     //    void cacheConvFunction(Int which, Array<Complex>& cf, CoordinateSystem& coord);
00396     // Find the convolution function
00397     void findConvFunction(const ImageInterface<Complex>& image,
00398                           const VisBuffer& vb);
00399     
00400     // Get the appropriate data pointer
00401     Array<Complex>* getDataPointer(const IPosition&, Bool);
00402     
00403     void ok();
00404     
00405     void init();
00406     //    virtual void initPolInfo(const VisBuffer& vb);
00407     // Is this record on Grid? check both ends. This assumes that the
00408     // ends bracket the middle
00409     Bool recordOnGrid(const VisBuffer& vb, Int rownr) const;
00410     
00411     // Padding in FFT
00412     Float padding_p;
00413     
00414     Int nWPlanes_p;
00415     // Image cache
00416     LatticeCache<Complex> * imageCache;
00417     
00418     // Sizes
00419     Long cachesize;
00420     Int tilesize;
00421     
00422     // Gridder
00423     ConvolveGridder<Double, Complex>* gridder;
00424     
00425     // Is this tiled?
00426     Bool isTiled;
00427     
00428     // Array lattice
00429     CountedPtr<Lattice<Complex> > arrayLattice;
00430     
00431     // Lattice. For non-tiled gridding, this will point to arrayLattice,
00432     //  whereas for tiled gridding, this points to the image
00433     CountedPtr<Lattice<Complex> > lattice;
00434     
00435     Float maxAbsData;
00436     
00437     // Useful IPositions
00438     IPosition centerLoc, offsetLoc;
00439     
00440     // Image Scaling and offset
00441     Vector<Double> uvScale, uvOffset;
00442     
00443     // Array for non-tiled gridding
00444     Array<Complex> griddedData;
00445     Array<DComplex> griddedData2;
00446     
00447     //    DirectionCoordinate directionCoord;
00448     MDirection::Convert* pointingToImage;
00449     
00450     // Grid/degrid zero spacing points?
00451     Bool usezero_p;
00452     
00453     //    CountedPtr<ConvolutionFunction> telescopeConvFunc_p;
00454     //    CFStore cfs_p, cfwts_p;
00455     // Array<Complex> convFunc_p, convWeights_p;
00456     //
00457     // Vector to hold the support size info. for the convolution
00458     // functions pointed to by the elements of convFunctions_p.  The
00459     // co-ordinates of this array are (W-term, Poln, PA).
00460     //
00461     Int convSize, convSampling, wConvSize, lastIndex_p;
00462     
00463     //
00464     // The average PB for sky image normalization
00465     //
00466     CountedPtr<ImageInterface<Float> > avgPB_p;
00467     CountedPtr<ImageInterface<Complex> > avgPBSq_p;
00468     //
00469     // No. of vis. polarization planes used in making the user defined
00470     // Stokes images
00471     //
00472     Int maxConvSupport;
00473     //
00474     // Percentage of the peak of the PB after which the image is set
00475     // to zero.
00476     //
00477     CountedPtr<SolvableVisJones> epJ_p;
00478     Double sigma;
00479     Int Nant_p, doPointing;
00480     Bool doPBCorrection, makingPSF, conjBeams_p;
00481     
00482     //    CountedPtr<CFCache> cfCache_p;
00483     ParAngleChangeDetector paChangeDetector;
00484     Double rotateOTFPAIncr_p, computePAIncr_p;
00485 
00486     Unit Second, Radian, Day;
00487     Array<Float> l_offsets,m_offsets;
00488     Vector<Float> pbPeaks, paList;
00489 
00490     Double currentCFPA, cfRefFreq_p, imRefFreq_p;
00491     Float lastPAUsedForWtImg;
00492     Bool pbNormalized_p;
00493     Vector<Bool> paNdxProcessed_p;
00494     //
00495     //----------------------------------------------------------------------
00496     //
00497     virtual void normalizeAvgPB();
00498     virtual void normalizeAvgPB(ImageInterface<Complex>& inImage,
00499                                 ImageInterface<Float>& outImage);
00500     virtual void resampleDataToGrid(Array<Complex>& griddedData, VBStore& vbs, 
00501                                     const VisBuffer& vb, Bool& dopsf);
00502     virtual void resampleDataToGrid(Array<DComplex>& griddedData, VBStore& vbs, 
00503                                     const VisBuffer& vb, Bool& dopsf);
00504     virtual void resampleGridToData(VBStore& vbs, Array<Complex>& griddedData,
00505                                     const VisBuffer& vb);
00506 
00507     virtual void setupVBStore(VBStore& vbs,
00508                               const VisBuffer& vb,
00509                               const Matrix<Float>& imagingweight,
00510                               const Cube<Complex>& visData,
00511                               const Matrix<Double>& uvw,
00512                               const Cube<Int>& flagCube,
00513                               const Vector<Double>& dphase);
00514 
00515     //    AWVisResampler visResampler_p;
00516     CountedPtr<VisibilityResamplerBase> visResampler_p;
00517     Int sensitivityPatternQualifier_p;
00518     String sensitivityPatternQualifierStr_p;
00519     CFStore rotatedConvFunc_p;
00520     CountedPtr<CFStore2> cfs2_p, cfwts2_p;
00521     Vector<Int> ConjCFMap_p, CFMap_p;
00522 
00523     Timer timer_p;
00524     Double runTime1_p;
00525 
00526 #include "AWProjectFT.FORTRANSTUFF.INC"
00527   };
00528 } //# NAMESPACE CASA - END
00529 
00530 #endif