FTMachine.h

Go to the documentation of this file.
00001 //# FTMachine.h: Definition for FTMachine
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_FTMACHINE_H
00030 #define SYNTHESIS_FTMACHINE_H
00031 
00032 #include <measures/Measures/Measure.h>
00033 #include <measures/Measures/MDirection.h>
00034 #include <measures/Measures/MPosition.h>
00035 #include <casa/Arrays/Array.h>
00036 #include <casa/Arrays/Vector.h>
00037 #include <casa/Arrays/Matrix.h>
00038 #include <casa/Logging/LogIO.h>
00039 #include <casa/Logging/LogSink.h>
00040 #include <casa/Logging/LogMessage.h>
00041 #include <casa/Containers/RecordInterface.h>
00042 #include <casa/Containers/Block.h>
00043 #include <images/Images/TempImage.h>
00044 #include <coordinates/Coordinates/SpectralCoordinate.h>
00045 #include <scimath/Mathematics/InterpolateArray1D.h>
00046 #include <synthesis/TransformMachines/CFCache.h>
00047 #include <synthesis/TransformMachines/CFStore2.h>
00048 
00049 #include <synthesis/TransformMachines/ConvolutionFunction.h>
00050 #include <synthesis/TransformMachines/PolOuterProduct.h>
00051 
00052 #include <images/Images/ImageInterface.h>
00053 #include <images/Images/SubImage.h>
00054 #include <synthesis/TransformMachines/StokesImageUtil.h>
00055 
00056 #include <synthesis/ImagerObjects/SIImageStore.h>
00057 #include <synthesis/ImagerObjects/SIImageStoreMultiTerm.h>
00058 
00059 namespace casa { //# NAMESPACE CASA - BEGIN
00060 
00061   class VisSet;
00062   class VisBuffer;
00063   class ROVisibilityIterator;
00064   class UVWMachine;
00065   class VisModelData;
00066   class SkyJones;
00067 
00068 // <summary> defines interface for the Fourier Transform Machine </summary>
00069 
00070 // <use visibility=export>
00071 
00072 // <reviewed reviewer="" date="" tests="" demos="">
00073 
00074 // <prerequisite>
00075 //   <li> <linkto class=SkyModel>SkyModel</linkto> module
00076 //   <li> <linkto class=SkyEquation>SkyEquation</linkto> module
00077 //   <li> <linkto class=VisBuffer>VisBuffer</linkto> module
00078 // </prerequisite>
00079 //
00080 // <etymology>
00081 // FTMachine is a Machine for Fourier Transforms
00082 // </etymology>
00083 //
00084 // <synopsis> 
00085 // The <linkto class=SkyEquation>SkyEquation</linkto> needs to be able
00086 // to perform Fourier transforms on visibility data. FTMachine
00087 // allows efficient Fourier Transform processing using a 
00088 // <linkto class=VisBuffer>VisBuffer</linkto> which encapsulates
00089 // a chunk of visibility (typically all baselines for one time)
00090 // together with all the information needed for processing
00091 // (e.g. UVW coordinates).
00092 // </synopsis> 
00093 //
00094 // <example>
00095 // A simple example of a FTMachine is found in 
00096 // <linkto class=GridFT>GridFT</linkto>.
00097 // See the example for <linkto class=SkyModel>SkyModel</linkto>.
00098 // </example>
00099 //
00100 // <motivation>
00101 // Define an interface to allow efficient processing of chunks of 
00102 // visibility data
00103 //
00104 // Note that the image must be Complex. It must contain the
00105 // Complex Stokes values (e.g. RR,RL,LR,LL). FTMachine
00106 // uses the image coordinate system to determine mappings
00107 // between the polarization and frequency values in the
00108 // PagedImage and in the VisBuffer.
00109 //
00110 // </motivation>
00111 //
00112 // <todo asof="97/10/01">
00113 // </todo>
00114 
00115 class FTMachine {
00116 public:
00117 
00118   //# Enumerations
00119   // Types of known Images that may be made using the makeImage method 
00120   enum Type {
00121     OBSERVED=0,         // From OBSERVED visibility data (default)
00122     MODEL,              // From MODEL visibility data
00123     CORRECTED,          // From CORRECTED visibility data
00124     RESIDUAL,           // From RESIDUAL (OBSERVED-MODEL) visibility data
00125     PSF,                // POINT SPREAD FUNCTION
00126     COVERAGE,           // COVERAGE (SD only)
00127     N_types,            // Number of types
00128     DEFAULT=OBSERVED
00129   };
00130 
00131   FTMachine();
00132 
00133 
00134   FTMachine(CountedPtr<CFCache>& cfcache,CountedPtr<ConvolutionFunction>& cfctor);
00135 
00136   FTMachine(const FTMachine& other);
00137 
00138   FTMachine& operator=(const FTMachine& other);
00139 
00140   void setBasePrivates(const FTMachine& other){FTMachine::operator=(other);}
00141 
00142   virtual ~FTMachine();
00143   
00144 
00145   //clone copy
00146   //should make it pure virtual forcing every ftm to have a cloner
00147   virtual FTMachine* cloneFTM(){return NULL;};
00148   // Initialize transform to Visibility plane
00149   virtual void initializeToVis(ImageInterface<Complex>& image, const VisBuffer& vb) = 0;
00150 
00151   // Vectorized InitializeToVis
00152   virtual void initializeToVis(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec,
00153                                PtrBlock<SubImage<Float> *> & modelImageVec, 
00154                                PtrBlock<SubImage<Float> *>& weightImageVec, 
00155                                PtrBlock<SubImage<Float> *>& fluxScaleVec, 
00156                                Block<Matrix<Float> >& weightsVec,
00157                                const VisBuffer& vb);
00158 
00159   virtual void initializeToVisNew(const VisBuffer& vb,
00160                                              CountedPtr<SIImageStore> imstore);
00161 
00162   //-------------------------------------------------------------------------------------
00163   // Finalize transform to Visibility plane
00164   // This is mostly a no-op, and is not-even called from CubeSkyEquation.
00165   virtual void finalizeToVis() = 0;
00166 
00167   // Note : No vectorized form of finalizeToVis yet.....
00168 
00169   //-------------------------------------------------------------------------------------
00170   // Initialize transform to Sky plane
00171   virtual void initializeToSky(ImageInterface<Complex>& image,
00172                                Matrix<Float>& weight, const VisBuffer& vb) = 0;
00173 
00174   // Vectorized InitializeToSky
00175   virtual void initializeToSky(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec,
00176                                Block<Matrix<Float> >& weightsVec, 
00177                                const VisBuffer& vb, 
00178                                const Bool dopsf=False);
00179 
00180   virtual void initializeToSkyNew(const Bool dopsf, 
00181                                   const VisBuffer& vb, 
00182                                   CountedPtr<SIImageStore> imstore);
00183 
00184   //-------------------------------------------------------------------------------------
00185   // Finalize transform to Sky plane
00186   virtual void finalizeToSky() = 0;
00187 
00188   virtual void finalizeToSky(ImageInterface<Complex>& iimage){(void)iimage;};
00189 
00190   // Vectorized finalizeToSky
00191   virtual void finalizeToSky(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec, 
00192                              PtrBlock<SubImage<Float> *> & resImageVec, 
00193                              PtrBlock<SubImage<Float> *>& weightImageVec, 
00194                              PtrBlock<SubImage<Float> *>& fluxScaleVec, 
00195                              Bool dopsf, 
00196                              Block<Matrix<Float> >& weightsVec, const VisBuffer& vb);
00197 
00198   virtual void finalizeToSkyNew(Bool dopsf, 
00199                                            const VisBuffer& vb,
00200                                            CountedPtr<SIImageStore> imstore  );
00201 
00202   //-------------------------------------------------------------------------------------
00203 
00204   // Get actual coherence from grid
00205   virtual void get(VisBuffer& vb, Int row=-1) = 0;
00206 
00207 
00208   // Put coherence to grid
00209   virtual void put(const VisBuffer& vb, Int row=-1, Bool dopsf=False, 
00210                    FTMachine::Type type= FTMachine::OBSERVED)=0;
00211 
00212   // Non const vb version - so that weights can be modified in-place
00213   // Currently, used only by MultiTermFT
00214   virtual void put(VisBuffer& vb, Int row=-1, Bool dopsf=False, 
00215                    FTMachine::Type type= FTMachine::OBSERVED)
00216                     {put((const VisBuffer&)vb,row,dopsf,type);};
00217 
00218   //-------------------------------------------------------------------------------------
00219   virtual void correlationToStokes(ImageInterface<Complex>& compImage, 
00220                                    ImageInterface<Float>& resImage, 
00221                                    const Bool dopsf);
00222  
00223   virtual void stokesToCorrelation(ImageInterface<Float>& modelImage,
00224                                    ImageInterface<Complex>& compImage);
00225 
00226   /*
00227   virtual void normalizeSumWeight(ImageInterface<Float>& inOutImage, 
00228                                ImageInterface<Float>& weightImage, 
00229                                const Bool dopsf);
00230   */
00231 
00232   virtual void normalizeImage(Lattice<Complex>&,//skyImage,
00233                               const Matrix<Double>&,// sumOfWts,
00234                               Lattice<Float>&,// sensitivityImage,
00235                               Bool /*fftNorm*/){return;};
00236 
00237   virtual void normalizeImage(ImageInterface<Float>& skyImage,
00238                               Matrix<Float>& sumOfWts,
00239                               ImageInterface<Float>& sensitivityImage,
00240                               Bool dopsf, Float pblimit, Int normtype);
00241 
00242 
00243   // All FTMachines that fill weightimage, need to set this.
00244   // TODO : Make this pure virtual.
00245   virtual Bool useWeightImage(){return False;}; 
00246   virtual Bool isSkyJonesSet(){return (sj_p.nelements()>0) && !(sj_p[0]).null()  ;}
00247   virtual Bool isSkyJonesChanged(VisBuffer& vb, Int row){if(sj_p.nelements()>0){return sj_p[0]->changed(vb,row);} else {return False;} };
00248 
00249   // Set SkyJones if image domain corrections /applycation are needed
00250   // To reset the the FTMachine for stopping image based correction/applycation
00251   // set in a Vector of size 0.
00252   // The pointers have to be handled by the caller ..no delete happening here
00253   virtual void setSkyJones(Vector<CountedPtr<SkyJones> >& sj);
00254   
00255   Bool changedSkyJonesLogic(const VisBuffer& vb, Bool& firstRow, Bool& internalRow);
00256 
00257 
00258   //-------------------------------------------------------------------------------------
00259 
00260   // Get the gridded visibilities or weight 
00261   template <typename T> void getGrid(Array<T>& thegrid);
00262    // Get the final image
00263   virtual ImageInterface<Complex>& getImage(Matrix<Float>&, Bool normalize=True) = 0;
00264 
00265   virtual const CountedPtr<ConvolutionFunction>& getAWConvFunc() {return convFuncCtor_p;};
00266 
00267   virtual void findConvFunction(const ImageInterface<Complex>&,// image,
00268                                 const VisBuffer& /*vb*/) {};
00269   // Get the final weights image
00270   virtual void getWeightImage(ImageInterface<Float>& weightImage, Matrix<Float>& weights) = 0;
00271 
00272   // Get a flux (divide by this to get a flux density correct image) 
00273   // image if there is one
00274   virtual void getFluxImage(ImageInterface<Float>& image){(void)image;};
00275 
00276   // Make the entire image
00277   virtual void makeImage(FTMachine::Type type,
00278                          VisSet& vs,
00279                          ImageInterface<Complex>& image,
00280                          Matrix<Float>& weight);
00281   // Make the entire image using a ROVisIter
00282   virtual void makeImage(FTMachine::Type type,
00283                          ROVisibilityIterator& vi,
00284                          ImageInterface<Complex>& image,
00285                          Matrix<Float>& weight);
00286 
00287   //-------------------------------------------------------------------------------------
00288 
00289   // Rotate the uvw from the observed phase center to the
00290   // desired phase center.
00291   void rotateUVW(Matrix<Double>& uvw, Vector<Double>& dphase,
00292                  const VisBuffer& vb);
00293 
00294   // Refocus on a finite distance
00295   void refocus(Matrix<Double>& uvw, const Vector<Int>& ant1,
00296                const Vector<Int>& ant2,
00297                Vector<Double>& dphase, const VisBuffer& vb);
00298 
00299   //helper function for openmp to call ...no private dependency
00300   static void locateuvw(const Double*& uvw, const Double*&dphase, const Double*& freq, const Int& nchan, const Double*& scale, const Double*& offset,  const Int& sampling, Int*& loc,Int*& off, Complex*& phasor, const Int& row, const Bool& doW=False); 
00301                  
00302 
00303   // Save and restore the FTMachine to and from a record
00304   virtual Bool toRecord(String& error, RecordInterface& outRecord, 
00305                         Bool withImage=False, const String diskimagename="");
00306   virtual Bool fromRecord(String& error, const RecordInterface& inRecord);
00307 
00308   // Has this operator changed since the last application?
00309   virtual Bool changed(const VisBuffer& vb);
00310 
00311   // Can this FTMachine be represented by Fourier convolutions?
00312   virtual Bool isFourier() {return False;}
00313 
00314   //set  spw for cube that will be used;
00315   Bool setSpw(Vector<Int>& spw, Bool validFrame);
00316 
00317   //return whether the ftmachine is using a double precision grid
00318   virtual Bool doublePrecGrid();
00319 
00320   // To make sure no padding is used in certain gridders
00321   virtual void setNoPadding(Bool nopad){(void)nopad;};
00322   
00323   // Return the name of the machine
00324 
00325   virtual String name() const =0;// { return "None";};
00326  
00327   // set and get the location used for frame 
00328   void setLocation(const MPosition& loc);
00329   MPosition& getLocation();
00330 
00331   // set a moving source aka planets or comets =>  adjust phase center
00332   // on the fly for gridding 
00333   virtual void setMovingSource(const String& sourcename);
00334   virtual void setMovingSource(const MDirection& mdir);
00335 
00336   //reset stuff in an FTMachine
00337   virtual void reset(){};
00338 
00339   //set frequency interpolation type
00340   virtual void setFreqInterpolation(const String& method);
00341 
00342   //tell ftmachine which Pointing table column to use for Direction
00343   //Mosaic or Single dish ft use this for example
00344   virtual void setPointingDirColumn(const String& column="DIRECTION");
00345 
00346   virtual String getPointingDirColumnInUse();
00347 
00348   virtual void setSpwChanSelection(const Cube<Int>& spwchansels);
00349   virtual Cube<Int> getSpwChanSelection(){return spwChanSelFlag_p;};
00350   virtual void setSpwFreqSelection(const Matrix<Double>& spwfreqs);
00351   virtual Matrix<Double> getSpwFreqSelection() {return spwFreqSel_p;};
00352 
00353   // set the order of the Taylor term for MFS this is to tell
00354   // A-Projection to qualify the accumulated avgPB for each Taylor
00355   // term in the CFCache.
00356   virtual void setMiscInfo(const Int qualifier)=0;
00357 
00358   virtual void setCanComputeResiduals(Bool& b) {canComputeResiduals_p=b;};
00359   virtual Bool canComputeResiduals() {return canComputeResiduals_p;};
00360   //
00361   // Make the VB and VBStore interefaces for the interim re-factoring
00362   // work.  Finally removed the VB interface.
00363   virtual void ComputeResiduals(VisBuffer&vb, Bool useCorrected) = 0;
00364   virtual Float getPBLimit() {return pbLimit_p;};
00365   //virtual void ComputeResiduals(VBStore& vb)=0;
00366   //get and set numthreads
00367   void setnumthreads(Int n);
00368   Int getnumthreads();
00369 
00370   virtual void setCFCache(CountedPtr<CFCache>& cfc, const Bool resetCFC=True);
00371   CountedPtr<CFCache> getCFCache() {return cfCache_p;};
00372   String getCacheDir() { return cfCache_p->getCacheDir(); };
00373 
00374   virtual void setDryRun(Bool val) 
00375   {
00376     isDryRun=val;
00377     //cerr << "FTM: " << isDryRun << endl;
00378   };
00379   virtual Bool dryRun() {return isDryRun;}
00380   virtual Bool isUsingCFCache() 
00381   {
00382     // cerr << "@#%$@% = " << cfCache_p.nrefs() << endl;
00383     return (cfCache_p.nrefs()!=0);
00384   }
00385   Bool isDryRun;
00386 
00387 protected:
00388 
00389   friend class VisModelData;
00390   friend class MultiTermFT;
00391   friend class MultiTermFTNew;
00392   LogIO logIO_p;
00393 
00394   LogIO& logIO();
00395 
00396   ImageInterface<Complex>* image;
00397 
00398   UVWMachine* uvwMachine_p;
00399 
00400   MeasFrame mFrame_p;
00401 
00402   // Direction of desired tangent plane
00403   Bool tangentSpecified_p;
00404   MDirection mTangent_p;
00405 
00406   MDirection mImage_p;
00407 
00408   // moving source stuff
00409   MDirection movingDir_p;
00410   Bool fixMovingSource_p;
00411   MDirection firstMovingDir_p;
00412     
00413 
00414   Double distance_p;
00415 
00416   uInt nAntenna_p;
00417 
00418   Int lastFieldId_p;
00419   Int lastMSId_p;
00420   //Use douple precision grid in gridding process
00421   Bool useDoubleGrid_p;
00422 
00423   virtual void initMaps(const VisBuffer& vb);
00424   virtual void initPolInfo(const VisBuffer& vb);
00425 
00426 
00427   // Sum of weights per polarization and per chan
00428   Matrix<Double> sumWeight, sumCFWeight;
00429 
00430   // Sizes
00431   Int nx, ny, npol, nchan, nvischan, nvispol;
00432 
00433   // Maps of channels and polarization
00434   Vector<Int> chanMap, polMap;
00435 
00436   // Is Stokes I only? iso XX,XY,YX,YY or LL,LR,RL,RR.
00437   Bool isIOnly;
00438 
00439   // Default Position used for phase rotations
00440   MPosition mLocation_p;
00441 
00442   // Set if uvwrotation is necessary
00443 
00444   Bool doUVWRotation_p;
00445   virtual void ok();
00446 
00447   // check if image is big enough for gridding
00448   
00449   virtual void gridOk (Int gridsupport);
00450 
00451   
00452   // setup multiple spectral window for cubes
00453   Block <Vector <Int> > multiChanMap_p;
00454   Vector<Int> selectedSpw_p;
00455   Vector<Int> nVisChan_p;
00456   Bool matchChannel(const Int& spw, 
00457                     const VisBuffer& vb);
00458 
00459   //redo all spw chan match especially if ms has changed underneath 
00460   Bool matchAllSpwChans(const VisBuffer& vb);
00461 
00462   //interpolate visibility data of vb to grid frequency definition
00463   //flag will be set the one as described in interpolateArray1D
00464   //return False if no interpolation is done...for e.g for nearest case
00465   virtual Bool interpolateFrequencyTogrid(const VisBuffer& vb, 
00466                                           const Matrix<Float>& wt,
00467                                           Cube<Complex>& data, 
00468                                           Cube<Int>& flag,
00469                                           Matrix<Float>& weight,
00470                                           FTMachine::Type type=FTMachine::OBSERVED ); 
00471   //degridded data interpolated back onto visibilities
00472   virtual Bool interpolateFrequencyFromgrid(VisBuffer& vb, 
00473                                             Cube<Complex>& data,
00474                                             FTMachine::Type type=FTMachine::MODEL );
00475   
00476 
00477   //Interpolate visibilities to be degridded upon
00478   virtual void getInterpolateArrays(const VisBuffer& vb,
00479                                     Cube<Complex>& data, Cube<Int>& flag);
00480 
00481 
00482   void setSpectralFlag(const VisBuffer& vb, Cube<Bool>& modflagcube); 
00483 
00484   //helper to save Measures in a record
00485   Bool saveMeasure(RecordInterface& rec, const String& name, String& error, const Measure& ms);
00486 
00487   // Private variables needed for spectral frame conversion 
00488   SpectralCoordinate spectralCoord_p;
00489   Vector<Bool> doConversion_p;
00490   Bool freqFrameValid_p;
00491   Vector<Double> imageFreq_p;
00492   //Vector of float lsrfreq needed for regridding
00493   Vector<Double> lsrFreq_p;
00494   Vector<Double> interpVisFreq_p;
00495   InterpolateArray1D<Double,Complex>::InterpolationMethod freqInterpMethod_p;
00496   String pointingDirCol_p;
00497   Cube<Int> spwChanSelFlag_p;
00498   Matrix<Double> spwFreqSel_p, expandedSpwFreqSel_p,expandedSpwConjFreqSel_p;
00499   Vector<Int> cfStokes_p;
00500   Int polInUse_p;
00501 
00502   CountedPtr<CFCache> cfCache_p;
00503   CFStore cfs_p, cfwts_p;
00504   CountedPtr<CFStore2> cfs2_p, cfwts2_p;
00505 
00506   CountedPtr<ConvolutionFunction> convFuncCtor_p;
00507   CountedPtr<PolOuterProduct> pop_p;
00508 
00509   Bool canComputeResiduals_p;
00510   Bool toVis_p;
00511   Int numthreads_p;
00512   
00513   // Array for non-tiled gridding
00514   // These are common to most FTmachines
00515   Array<Complex> griddedData;
00516   Array<DComplex> griddedData2;
00517 
00518 
00519   Float pbLimit_p;
00520   //  Vector<SkyJones *> sj_p;
00521   Vector<CountedPtr<SkyJones> > sj_p;
00522   //A holder for the complex image if nobody else is keeping it
00523   CountedPtr<ImageInterface<Complex> > cmplxImage_p;
00524 
00525  private:
00526   //Some temporary wasteful function for swapping axes because we don't 
00527   //Interpolation along the second axis...will need to implement 
00528   //interpolation on y axis of a cube. 
00529   
00530   void swapyz(Cube<Complex>& out, const Cube<Complex>& in);
00531   void swapyz(Cube<Complex>& out, const Cube<Bool>& outFlag, const Cube<Complex>& in);
00532   void swapyz(Cube<Bool>& out, const Cube<Bool>& in);
00533   void convUVW(Double& dphase, Vector<Double>& thisrow);
00534 
00535 };
00536 
00537 #include <synthesis/TransformMachines/FTMachine.tcc>
00538 
00539 
00540 } //# NAMESPACE CASA - END
00541 
00542 #endif
00543 
00544 
00545 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1