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

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1