CFBuffer.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //# CFBuffer.h: Definition of the CFBuffer class
00003 //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
00004 //# Associated Universities, Inc. Washington DC, USA.
00005 //#
00006 //# This library is free software; you can redistribute it and/or modify it
00007 //# under the terms of the GNU Library General Public License as published by
00008 //# the Free Software Foundation; either version 2 of the License, or (at your
00009 //# option) any later version.
00010 //#
00011 //# This library is distributed in the hope that it will be useful, but WITHOUT
00012 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00013 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00014 //# License for more details.
00015 //#
00016 //# You should have received a copy of the GNU Library General Public License
00017 //# along with this library; if not, write to the Free Software Foundation,
00018 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00019 //#
00020 //# Correspondence concerning AIPS++ should be addressed as follows:
00021 //#        Internet email: aips2-request@nrao.edu.
00022 //#        Postal address: AIPS++ Project Office
00023 //#                        National Radio Astronomy Observatory
00024 //#                        520 Edgemont Road
00025 //#                        Charlottesville, VA 22903-2475 USA
00026 //#
00027 //# $Id$
00028 #ifndef SYNTHESIS_CFBUFFER_H
00029 #define SYNTHESIS_CFBUFFER_H
00030 #include <synthesis/TransformMachines/SynthesisError.h>
00031 #include <coordinates/Coordinates/CoordinateSystem.h>
00032 #include <synthesis/TransformMachines/CFDefs.h>
00033 #include <synthesis/TransformMachines/CFCell.h>
00034 #include <synthesis/TransformMachines/Utils.h>
00035 #include <images/Images/ImageInterface.h>
00036 #include <casa/Utilities/CountedPtr.h>
00037 #include <casa/Utilities/Sort.h>
00038 #include <casa/Logging/LogOrigin.h>
00039 #include <msvis/MSVis/VisBuffer.h>
00040 #include <casa/Logging/LogSink.h>
00041 #include <casa/Logging/LogIO.h>
00042 //
00043 // <summary> defines interface for the storage for convolution functions </summary>
00044 // <use visibility=export>
00045 //
00046 // <prerequisite>
00047 // </prerequisite>
00048 //
00049 // <etymology> 
00050 //
00051 //  CFBuffer is the basic in-memory storage for convolution functions
00052 //  as a function of polarization, W-value and frequency at a
00053 //  particular value of Parallactic Angle and baseline type.
00054 //
00055 //</etymology>
00056 //
00057 // <synopsis> 
00058 //
00059 // The CFBuffer class encapsulates the storage and associated
00060 // auxillary information required for the convolution functions.  The
00061 // <linto class=CFStore>CFStore</linkto> class is a collection of
00062 // CFBuffer objects.  A collection of CFStore objects is held and
00063 // managed by the <linto class=FTMachine>FTMachine</linkto> and the
00064 // appropriate one, depending on the Time/PA value, polarization and
00065 // frequency of the data in the <linkto
00066 // class=VisBuffer>VisBuffer</linkto>, is supplied to the <linkto
00067 // class=VisibilityResampler>VisibilityResampler</linkto> object for
00068 // re-sampling the data onto a grid (or vice versa).
00069 //
00070 // Conceptually, this object holds the convolution functions
00071 // parameterized by the properties of the electromagnetic radiation
00072 // (polarization state, frequency and the w-value which is the Fresnel
00073 // term and implicitly a function of the frequency).  The <linkto
00074 // class=CFStore>CFStore</linkto> object holds a list of this object
00075 // index by the telescope related parameters (antenna1,
00076 // antenna2, Parallactic Angle or Time, etc.).
00077 //
00078 //</synopsis>
00079 //
00080 // <example>
00081 // </example>
00082 //
00083 // <motivation>
00084 //
00085 // To factor out the details of effecient in-memory storage of the
00086 // convolution functions into a separate class.  This class can then
00087 // be optmized by specializations for various types of convolution
00088 // functions and imaging algorithms without the need to change the
00089 // imaging framework.
00090 //
00091 // </motivation>
00092 //
00093 
00094 using namespace casa::CFDefs;
00095 using namespace std;
00096 namespace casa { //# NAMESPACE CASA - BEGIN
00097   //  template <class T>
00098   typedef Complex TT;
00099 
00100   struct CFBStruct {
00101     CFCStruct *CFBStorage;
00102     int shape[3];
00103     Double *freqValues, *wValues, *pointingOffset;
00104     Double fIncr, wIncr;
00105     Int **muellerElementsIndex, **muellerElements, 
00106       **conjMuellerElementsIndex, **conjMuellerElements,
00107       **conjFreqNdxMap, **freqNdxMap;
00108     Int nMueller;
00109     
00110 
00111     CFCStruct* getCFB(int i, int j, int k)
00112     { return &(CFBStorage[i + (shape[1]-1)*j + (shape[1]-1)*(shape[2]-1)*k]);}
00113     //    { return &(CFBStorage[i + (shape[1]-1)*j + (shape[2]-1)*k]);}
00114   } ;
00115 
00116   class CFBuffer
00117   {
00118   public:
00119     //
00120     //========================= Administrative Parts ==========================
00121     //------------------------------------------------------------------
00122     //
00123     CFBuffer(): wValues_p(), maxXSupport_p(-1), maxYSupport_p(-1), pointingOffset_p(), cfHitsStats(),
00124                 freqNdxMapsReady_p(False), freqNdxMap_p(), conjFreqNdxMap_p()
00125     {};
00126     
00127     CFBuffer(Int maxXSup, Int maxYSup):
00128       wValues_p(), maxXSupport_p(maxXSup), maxYSupport_p(maxYSup), pointingOffset_p(), cfHitsStats(),
00129       freqNdxMapsReady_p(False), freqNdxMap_p(), conjFreqNdxMap_p()
00130     {
00131       // storage_p.resize(1,1,1); 
00132       // storage_p(0,0,0) = new Array<TT>(dataPtr);
00133       // coordSys_p.resize(1,1,1); 
00134       // coordSys_p(0,0,0) = cs;
00135     };
00136     
00137     ~CFBuffer() 
00138     {
00139       //cerr << "############### " << "~CFBuffer() called" << endl;
00140       // LogIO log_l(LogOrigin("CFBuffer","~CFBuffer[R&D]"));
00141       // log_l << "CF Hits stats gathered: " << cfHitsStats << endl;
00142     };
00143     
00144     CountedPtr<CFBuffer> clone();
00145     void allocCells(const Cube<CountedPtr<CFCell> >& cells);
00146     void setParams(const CFBuffer& other);
00147     //
00148     //============================= Functional Parts ============================
00149     //------------------------------------------------------------------
00150     //
00151     //    CFBuffer& operator=(const CFBuffer& other);
00152     //
00153     // Get the single convolution function as an Array<T> for the
00154     // supplied value of the frequency and the muellerElement.
00155     // Mueller element is essentially the polarization product, but
00156     // can be any of the of 16 elements of the outer product.
00157     //
00158     //-------------------------------------------------------------------------
00159     //
00160     inline Int nChan() {return nChan_p;}
00161     inline Int nW() {return nW_p;}
00162     inline Int nMuellerElements() {return nPol_p;}
00163     inline IPosition shape() {IPosition shp(3,nChan_p, nW_p, nPol_p); return shp;}
00164     
00165     inline Vector<Double> getFreqList() {return freqValues_p;};
00166     inline Vector<Double> getWList() {return wValues_p;};
00167     
00168     CFCell& getCFCell(const Double& freqVal, const Double& wValue, 
00169                       const Int & muellerElement); 
00170     // muellerElement: (i,j) of the Mueller Matrix
00171     CountedPtr<CFCell>& getCFCellPtr(const Double& freqVal, const Double& wValue, 
00172                                      const Int & muellerElement); 
00173     CFCell& operator()(const Int& i, const Int& j, const Int& k) {return *cfCells_p(i,j,k);}
00174     CFCell& getCFCell(const Int& i, const Int& j, const Int& k);
00175 
00176     CountedPtr<CFCell >& getCFCellPtr(const Int& i, const Int& j, const Int& k);
00177     
00178     //=========================================================================
00179     Array<TT>& getCF(const Double& freqVal, const Double& wValue, 
00180                      const Int & muellerElement)
00181     {return *(getCFCell(freqVal, wValue, muellerElement).storage_p);}
00182     // muellerElement: (i,j) of the Mueller Matrix
00183     
00184     CountedPtr<Array<TT> >& getCFPtr(const Double& freqVal, const Double& wValue, 
00185                                      const Int & muellerElement) 
00186     {return getCFCellPtr(freqVal, wValue, muellerElement)->storage_p;}
00187     
00188     Array<TT>& getCF(const Int& i, const Int& j, const Int& k)
00189     {return *(getCFCell(i,j,k).storage_p);}
00190     
00191     CountedPtr<Array<TT> >& getCFPtr(const Int& i, const Int& j, const Int& k)
00192     {return getCFCellPtr(i,j,k)->storage_p;}
00193     
00194     
00195     //
00196     // Get the parameters of a the CFs indexed by values.  The version
00197     // which returns also the Coordinate System associated with the
00198     // CFs are slow (CoordinateSystem::operator=() is surprisingly
00199     // expensive!).  So do not use this in tight loops.  If it is
00200     // required, use the version without the co-ordinate system below.
00201     //
00202     void getParams(CoordinateSystem& cs, Float& sampling, 
00203                    Int& xSupport, Int& ySupport, 
00204                    const Double& freqVal, const Double& wValue, 
00205                    const Int& muellerElement);
00206     //-------------------------------------------------------------------------
00207     // Get CF by directly indexing in the list of CFs (data vector)
00208     inline void getParams(CoordinateSystem& cs, Float& sampling, 
00209                           Int& xSupport, Int& ySupport, 
00210                           const Int& i, const Int& j, const Int& k)
00211     {
00212       cs = cfCells_p(i,j,k)->coordSys_p;
00213       sampling = cfCells_p(i,j,k)->sampling_p;
00214       xSupport = cfCells_p(i,j,k)->xSupport_p;
00215       ySupport = cfCells_p(i,j,k)->ySupport_p;
00216     }
00217     void getParams(Double& freqVal, Float& sampling, 
00218                    Int& xSupport, Int& ySupport, 
00219                    const Int& iFreq, const Int& iW, const Int& iPol)
00220     {
00221       sampling = cfCells_p(iFreq,iW,iPol)->sampling_p;
00222       xSupport = cfCells_p(iFreq,iW,iPol)->xSupport_p;
00223       ySupport = cfCells_p(iFreq,iW,iPol)->ySupport_p;
00224       freqVal = freqValues_p(iFreq);
00225     }
00226     
00227     inline void getCoordList(Vector<Double>& freqValues, Vector<Double>& wValues,
00228                              PolMapType& muellerElementsIndex, PolMapType& muellerElements, 
00229                              PolMapType& conjMuellerElementsIndex, PolMapType& conjMuellerElements, 
00230                              Double& fIncr, Double& wIncr)
00231     {
00232       freqValues.assign(freqValues_p);wValues.assign(wValues_p);
00233       muellerElements.assign(muellerElements_p);         muellerElementsIndex.assign(muellerElementsIndex_p);
00234       conjMuellerElements.assign(conjMuellerElements_p); conjMuellerElementsIndex.assign(conjMuellerElementsIndex_p);
00235       fIncr = freqValIncr_p; wIncr = wValIncr_p;
00236     }
00237     
00238     Int nearestNdx(const Double& val, const Vector<Double>& valList, const Double& incr);
00239     
00240     Int nearestFreqNdx(const Double& freqVal) ;
00241     
00242     inline Int nearestWNdx(const Double& wVal) 
00243     {
00244       //      return SynthesisUtils::nint(sqrt(wValIncr_p*abs(wVal)));
00245       return max(0,min((int)(sqrt(wValIncr_p*abs(wVal))),wValues_p.nelements())-1);
00246       // Int ndx=(int)(sqrt(wValIncr_p*abs(wVal)));
00247       // if ((uInt)ndx >= wValues_p.nelements())
00248       //        cerr << endl << endl << ndx << " " <<  wVal << " " << wValIncr_p << endl << endl;
00249       // return min(ndx,wValues_p.nelements()-1);
00250     }
00251     
00252     Double nearest(Bool& found, const Double& val, const Vector<Double>& valList, const Double& incr);
00253     
00254     inline Double nearestFreq(Bool& found, const Double& freqVal)
00255     {return nearest(found, freqVal, freqValues_p, freqValIncr_p);}
00256     
00257     inline Double nearestWVal(Bool& found, const Double& wVal)
00258     {return nearest(found, wVal, wValues_p, wValIncr_p);}
00259     
00260     //-------------------------------------------------------------------------
00261     //
00262     // Generate a map for the given frequency and Mueller element list
00263     // to the index in the internal list of CFs.  This can be used in
00264     // tight loops to get get direct access to the required CF.
00265     //
00266     void makeCFBufferMap(const Vector<Double>& freqVals, 
00267                          const Vector<Double>& wValues,
00268                          const MuellerMatrixType& muellerElements);
00269     //-------------------------------------------------------------------------
00270     //
00271     // Add a Convolution Function with associated parameters.
00272     //
00273     void addCF(Array<TT>*, //dataPtr, 
00274                CoordinateSystem&,// cs, 
00275                Float& ,//sampling, 
00276                Int& ,//xSupport, 
00277                Int& ,//ySupport,
00278                Double& ,//freqValue, 
00279                Double& ,//wValue, 
00280                Int& //muellerElement
00281                )
00282     {throw(AipsError("CFBuffer::addCF called"));}
00283     //-------------------------------------------------------------------------
00284     //
00285     void resize(const IPosition& size) {cfCells_p.resize(size);};
00286     void resize(const Double& wIncr, const Double& freqIncr,
00287                 const Vector<Double>& wValues, 
00288                 const Vector<Double>& freqValues,
00289                 const PolMapType& muellerElements,
00290                 const PolMapType& muellerElementsIndex,
00291                 const PolMapType& conjMuellerElements,
00292                 const PolMapType& conjMuellerElementsIndex);
00293     Int noOfMuellerElements(const PolMapType& muellerElements);
00294     //-------------------------------------------------------------------------
00295     // Set only the CF parameters.  Return to index of the CF that was set.
00296     //
00297     RigidVector<Int, 3> setParams(const Int& i, const Int& j, const Int& ipx, const Int& ipy,
00298                                   CoordinateSystem& cs, Float& sampling,
00299                                   Int& xSupport, Int& ySupport,
00300                                   const Double& freqValue, const Double& wValue, 
00301                                   const Int& muellerElement,
00302                                   const String& fileName=String(),
00303                                   const Double& conjFreq=0.0,
00304                                   const Int& conjPol=-1,
00305                                   const String& telescopeName=String(),
00306                                   const Float& diameter=25.0);
00307     // RigidVector<Int, 3> setParams(const Int& inu, const Int& iw, const Int& muellerElement,
00308     //                            const TableRecord& miscInfo);
00309     void setPointingOffset(const Vector<Double>& offset) 
00310     {pointingOffset_p.assign(offset);};
00311     Vector<Double> getPointingOffset() {return pointingOffset_p;};
00312     //
00313     // Also set the size of the CF in x and y.
00314     //
00315     void setParams(Int& nx, Int& ny, CoordinateSystem& cs, Float& sampling, 
00316                    Int& xSupport, Int& ySupport, 
00317                    const Double& freqVal, const Double& wValue, 
00318                    const Int& muellerElement,
00319                    const String& fileName);
00320     void setPA(Float& pa);
00321     RigidVector<Int,3> getIndex(const Double& freqVal, const Double& wValue, 
00322                                 const Int& muellerElement);
00323     //-------------------------------------------------------------------------
00324     //
00325     // Copy just the parameters from other to this.
00326     //
00327     void copyParams(const CFBuffer& other)
00328     {
00329       cfCells_p = other.cfCells_p;
00330       // coordSys_p = other.coordSys_p; sampling_p.assign(other.sampling_p); 
00331       // xSupport_p.assign(other.xSupport_p); ySupport_p.assign(other.ySupport_p);
00332       maxXSupport_p=other.maxXSupport_p;  maxYSupport_p=other.maxYSupport_p; 
00333     }
00334     //-------------------------------------------------------------------------
00335     //
00336     // Write the description of the storage on the supplied ostream.
00337     // Used mostly for debugging, but might be useful for user
00338     // feedback/logging.
00339     //
00340     void show(const char *Mesg=NULL,ostream &os=cerr);
00341     //
00342     // Returns True if the internal storage is not yet initialized.
00343     //
00344     Bool null() {return (cfCells_p.nelements() == 0);};
00345     
00346     Cube<CountedPtr<CFCell> >& getStorage() {return cfCells_p;};
00347     void makePersistent(const char *dir, const char *cfName="");
00348     
00349     void primeTheCache();
00350     void initMaps(const VisBuffer& vb,const Matrix<Double>& freqSelection,const Double& imRefFreq);
00351     void initPolMaps(PolMapType& polMap, PolMapType& conjPolMap);
00352     //
00353     // For CUDA kernel
00354     //
00355     void getFreqNdxMaps(Vector<Vector<Int> >& freqNdx, Vector<Vector<Int> >& conjFreqNdx);
00356     inline Int nearestFreqNdx(const Int& spw, const Int& chan, const Bool conj=False)
00357     {
00358       if (conj) return conjFreqNdxMap_p[spw][chan];
00359       else  return freqNdxMap_p[spw][chan];
00360     }
00361     
00362     void getAsStruct(CFBStruct& st);
00363     
00364     static void initCFBStruct(CFBStruct& cfbSt) 
00365     {
00366       cfbSt.CFBStorage=NULL;
00367       cfbSt.freqValues=NULL;
00368       cfbSt.wValues=NULL;
00369       cfbSt.muellerElementsIndex=NULL;
00370       cfbSt.muellerElements=NULL;
00371       cfbSt.conjMuellerElementsIndex=NULL;
00372       cfbSt.conjMuellerElements=NULL;
00373       cfbSt.shape[0]=cfbSt.shape[1]=cfbSt.shape[2]=0;
00374       cfbSt.fIncr=cfbSt.wIncr=0.0;
00375     }
00376     void fill(const Int& nx, const Int& ny, 
00377               const Vector<Double>& freqValues,
00378               const Vector<Double>& wValues,
00379               const PolMapType& muellerElements);
00380     
00381     IPosition getShape() {return cfCells_p.shape();}
00382     //
00383     //============================= Protected Parts ============================
00384     //------------------------------------------------------------------
00385     //
00386   protected:
00387     //
00388     // The storage buffer for the pixel values in CFCell is Array<T>
00389     // rather than Matrix<T> to accomodate rotationally symmetric CFs
00390     // (like the Prolate Spheroidal) which can be held as a Vector of
00391     // values.
00392     //
00393     Cube<CountedPtr<CFCell> > cfCells_p;// freqValues x wValues x muellerElements
00394     Vector<Double> wValues_p, freqValues_p;
00395     PolMapType muellerElements_p, muellerElementsIndex_p,conjMuellerElements_p,conjMuellerElementsIndex_p; 
00396     Double wValIncr_p, freqValIncr_p;
00397     MuellerMatrixType muellerMask_p;
00398     
00399     Int nPol_p, nChan_p, nW_p, maxXSupport_p, maxYSupport_p;
00400     Vector<Double> pointingOffset_p;
00401     Cube<Int> cfHitsStats;
00402     Bool freqNdxMapsReady_p;
00403     Vector<Vector<Int> > freqNdxMap_p, conjFreqNdxMap_p;
00404     void ASSIGNVVofI(Int** &target,Vector<Vector<Int> >& source, Bool& doAlloc);
00405   };
00406 
00407   // declare a commonly used template extern
00408   extern template class casa::Array<casa::CountedPtr<casa::CFBuffer> >;
00409   
00410 } //# NAMESPACE CASA - END
00411 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1