AWVisResampler.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //# AWVisResampler.h: Definition of the AWVisResampler 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 
00029 #ifndef SYNTHESIS_AWVISRESAMPLER_H
00030 #define SYNTHESIS_AWVISRESAMPLER_H
00031 
00032 #include <synthesis/TransformMachines/CFStore.h>
00033 #include <synthesis/TransformMachines/VBStore.h>
00034 #include <synthesis/TransformMachines/VisibilityResampler.h>
00035 #include <msvis/MSVis/VisBuffer.h>
00036 #include <casa/Arrays/Array.h>
00037 #include <casa/Arrays/Vector.h>
00038 
00039 #include <casa/Logging/LogIO.h>
00040 #include <casa/Logging/LogSink.h>
00041 #include <casa/Logging/LogMessage.h>
00042 
00043 namespace casa { //# NAMESPACE CASA - BEGIN
00044   class AWVisResampler: public VisibilityResampler
00045   {
00046   public: 
00047     AWVisResampler(): VisibilityResampler(),
00048                       cached_phaseGrad_p(),
00049                       cached_PointingOffset_p()
00050     {cached_PointingOffset_p.resize(2);cached_PointingOffset_p=-1000.0;runTimeG_p=runTimeDG_p=0.0;};
00051     //    AWVisResampler(const CFStore& cfs): VisibilityResampler(cfs)      {}
00052     virtual ~AWVisResampler()                                         {};
00053 
00054     virtual VisibilityResamplerBase* clone()
00055     {return new AWVisResampler(*this);}
00056     
00057     // AWVisResampler(const AWVisResampler& other): VisibilityResampler(other),cfMap_p(), conjCFMap_p()
00058     // {copy(other);}
00059 
00060     virtual void copyMaps(const AWVisResampler& other)
00061     {setCFMaps(other.cfMap_p, other.conjCFMap_p);}
00062     virtual void copy(const VisibilityResamplerBase& other) 
00063     {
00064       VisibilityResampler::copy(other);
00065       // const Vector<Int> cfmap=other.getCFMap();
00066       // const Vector<Int> conjcfmap = other.getConjCFMap();
00067 
00068       // setCFMaps(cfmap,conjcfmap);
00069     }
00070 
00071     virtual void copy(const AWVisResampler& other) 
00072     {
00073       VisibilityResampler::copy(other);
00074       SynthesisUtils::SETVEC(cached_phaseGrad_p, other.cached_phaseGrad_p);
00075       SynthesisUtils::SETVEC(cached_PointingOffset_p, other.cached_PointingOffset_p);
00076     }
00077 
00078     AWVisResampler& operator=(const AWVisResampler& other) 
00079     {
00080       copy(other);      
00081       SynthesisUtils::SETVEC(cached_phaseGrad_p, other.cached_phaseGrad_p);
00082       SynthesisUtils::SETVEC(cached_PointingOffset_p, other.cached_PointingOffset_p);
00083       return *this;
00084     }
00085 
00086     virtual void setCFMaps(const Vector<Int>& cfMap, const Vector<Int>& conjCFMap)
00087     {SETVEC(cfMap_p,cfMap);SETVEC(conjCFMap_p,conjCFMap);}
00088 
00089     // virtual void setConvFunc(const CFStore& cfs) {convFuncStore_p = cfs;};
00090     //
00091     //------------------------------------------------------------------------------
00092     //
00093     // Re-sample the griddedData on the VisBuffer (a.k.a gridding).
00094     //
00095     // In this class, these just call the private templated version.
00096     // The first variant grids onto a double precision grid while the
00097     // second one does it on a single precision grid.
00098     //
00099     // Note that the following calls allow using any CFStore object
00100     // for gridding while de-gridding uses the internal
00101     // convFuncStore_p object.
00102     // virtual void DataToGrid(Array<DComplex>& griddedData, VBStore& vbs, Matrix<Double>& sumwt,
00103     //                      const Bool& dopsf, CFStore& cfs)
00104     // {DataToGridImpl_p(griddedData, vbs, sumwt,dopsf,cfs);}
00105 
00106     // virtual void DataToGrid(Array<Complex>& griddedData, VBStore& vbs, Matrix<Double>& sumwt,
00107     //                      const Bool& dopsf, CFStore& cfs)
00108     // {DataToGridImpl_p(griddedData, vbs, sumwt,dopsf,cfs);}
00109     //
00110     // Simulating defaulting CFStore arguemnt in the above calls to convFuncStore_p
00111     //
00112 
00113     //***TEMP REMOVAL OF DComplex gridder*****
00114 
00115     virtual void DataToGrid(Array<DComplex>& griddedData, VBStore& vbs, Matrix<Double>& sumwt,
00116                             const Bool& dopsf,Bool useConjFreqCF=False)
00117     {DataToGridImpl_p(griddedData, vbs, sumwt,dopsf,useConjFreqCF);}
00118 
00119     virtual void DataToGrid(Array<Complex>& griddedData, VBStore& vbs, Matrix<Double>& sumwt,
00120                             const Bool& dopsf,Bool useConjFreqCF=False)
00121     {DataToGridImpl_p(griddedData, vbs, sumwt,dopsf,useConjFreqCF);}
00122 
00123     //
00124     //------------------------------------------------------------------------------
00125     //
00126     // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding)
00127     //
00128     virtual void GridToData(VBStore& vbs,const Array<Complex>& griddedData); 
00129     //    virtual void GridToData(VBStore& vbs, Array<Complex>& griddedData); 
00130   protected:
00131     virtual Complex getConvFuncVal(const Cube<Double>& convFunc, const Matrix<Double>& uvw, 
00132                                    const Int& irow, const Vector<Int>& pixel)
00133     {
00134       (void)uvw; (void)irow;return convFunc(pixel[0],pixel[1],pixel[2]);
00135     }
00136     Complex getCFArea(Complex* __restrict__& convFuncV, Double& wVal,
00137                       Vector<Int>& scaledSupport, Vector<Float>& scaledSampling,
00138                       Vector<Double>& off,
00139                       Vector<Int>& convOrigin, Vector<Int>& cfShape,
00140                       Double& sinDPA, Double& cosDPA);
00141 
00142   template <class T>
00143   Complex accumulateOnGrid(Array<T>& grid, Complex* __restrict__& convFuncV, 
00144                            Complex& nvalue,
00145                            Double& wVal, Vector<Int>& scaledSupport, 
00146                            Vector<Float>& scaledSampling, Vector<Double>& off,
00147                            Vector<Int>& convOrigin, Vector<Int>& /*cfShape*/,
00148                            Vector<Int>& loc, Vector<Int>& igrdpos, 
00149                            Double& /*sinDPA*/, Double& /*cosDPA*/,
00150                            Bool& finitePointingOffset, Bool dopsf);
00151   template <class T>
00152   void XInnerLoop(const Int *scaleSupport, const Float* scaledSampling,
00153                   const Double* off,
00154                   const Int* loc, Complex& cfArea,  
00155                   const Int * __restrict__ iGrdPosPtr,
00156                   Complex *__restrict__& convFuncV,
00157                   const Int* convOrigin,
00158                   Complex& nvalue,
00159                   Double& wVal,
00160                   Bool& /*finitePointingOffset*/,
00161                   Bool& /*doPSFOnly*/,
00162                   T* __restrict__ gridStore,
00163                   Int* iloc,
00164                   Complex& norm,
00165                   Int* igrdpos);
00166 
00167   template <class T>
00168   void accumulateFromGrid(T& nvalue, const T* __restrict__& grid, 
00169                           Vector<Int>& iGrdPos,
00170                           Complex* __restrict__& convFuncV, 
00171                           Double& wVal, Vector<Int>& scaledSupport, 
00172                           Vector<Float>& scaledSampling, Vector<Double>& off,
00173                           Vector<Int>& convOrigin, Vector<Int>& cfShape,
00174                           Vector<Int>& loc, 
00175                           Complex& phasor, 
00176                           Double& sinDPA, Double& cosDPA,
00177                           Bool& finitePointingOffset, 
00178                           Matrix<Complex>& cached_phaseGrad_p);
00179 
00180     //
00181     //------------------------------------------------------------------------------
00182     //----------------------------Private parts-------------------------------------
00183     //------------------------------------------------------------------------------
00184     //
00185   private:
00186     // Vector<Double> uvwScale_p, offset_p, dphase_p;
00187     // Vector<Int> chanMap_p, polMap_p;
00188     // CFStore convFuncStore_p;
00189     // //    Int inc0_p, inc1_p, inc2_p, inc3_p;
00190     // Vector<Int> inc_p;
00191     //    Vector<Int> cfMap_p, conjCFMap_p;
00192     Vector<Int> gridInc_p, cfInc_p;
00193     Matrix<Complex> cached_phaseGrad_p;
00194     Vector<Double> cached_PointingOffset_p;
00195     //
00196     // Re-sample the griddedData on the VisBuffer (a.k.a de-gridding).
00197     //
00198     template <class T>
00199     void DataToGridImpl_p(Array<T>& griddedData, VBStore& vb,  
00200                           Matrix<Double>& sumwt,const Bool& dopsf,
00201                           Bool /*useConjFreqCF*/);
00202 
00203     void sgrid(Vector<Double>& pos, Vector<Int>& loc, Vector<Double>& off, 
00204                Complex& phasor, const Int& irow, const Matrix<Double>& uvw, 
00205                const Double& dphase, const Double& freq, 
00206                const Vector<Double>& scale, const Vector<Double>& offset,
00207                const Vector<Float>& sampling);
00208 
00209     inline Bool onGrid (const Int& nx, const Int& ny, const Int& nw, 
00210                         const Vector<Int>& loc, 
00211                         const Vector<Int>& support)
00212     {
00213       return (((loc(0)-support[0]) >= 0 ) && ((loc(0)+support[0]) < nx) &&
00214               ((loc(1)-support[1]) >= 0 ) && ((loc(1)+support[1]) < ny) &&
00215               (loc(2) >= 0) && (loc(2) <= nw));
00216     };
00217 
00218     // Array assignment operator in CASACore requires lhs.nelements()
00219     // == 0 or lhs.nelements()=rhs.nelements()
00220     template <class T>
00221     inline void SETVEC(Vector<T>& lhs, const Vector<T>& rhs)
00222     {lhs.resize(rhs.shape()); lhs = rhs;};
00223 
00224 
00225     //
00226     // Internal methods to address a 4D array.  These should ulimately
00227     // moved to a Array4D class in CASACore
00228     //
00229 
00230     // This is called less frequently.  Currently once per VisBuffer
00231     // inline void cacheAxisIncrements(const Vector<Int>& n, Vector<Int>& inc)
00232     // {inc.resize(4);inc[0]=1, inc[1]=inc[0]*n[0], inc[2]=inc[1]*n[1], inc[3]=inc[2]*n[2];(void)n[3];}
00233 
00234 
00235     // The following method is also called from the inner loop, but
00236     // does not use CASA Vector (which are not thread safe, I (SB) am
00237     // told).
00238     inline Complex getFrom4DArray(const Complex *__restrict__& store,
00239                                   const Int* iPos, const Int* inc)
00240     {
00241       return *(store+(iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]));
00242       //      return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];
00243     };
00244 
00245     // The following two methods are called in the innermost loop.
00246     inline Complex getFrom4DArray(const Complex *__restrict__& store,
00247                                   const Vector<Int>& iPos, const Vector<Int>& inc)
00248     {
00249       return *(store+(iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]));
00250       //      return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];
00251     };
00252     inline DComplex getFrom4DArray(const DComplex *__restrict__& store,
00253                                   const Vector<Int>& iPos, const Vector<Int>& inc)
00254     {
00255       return *(store+(iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]));
00256       //      return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];
00257     };
00258 
00259     template <class T>
00260     void addTo4DArray(T *__restrict__& store,
00261                       const Int *__restrict__& iPos, const Vector<Int>& inc, 
00262                       Complex& nvalue, Complex& wt) __restrict__
00263     {
00264       // T *tmp=store+(iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]);
00265       // *tmp += nvalue*wt;
00266       store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]] += (nvalue*wt);
00267     }
00268 
00269     //
00270     // This rotates the convolution function by rotating the
00271     // co-ordinate system.  For the accuracies already required for
00272     // EVLA and ALMA, this is not useful.  Leaving it hear for now....
00273     //
00274     Bool reindex(const Vector<Int>& in, Vector<Int>& out,
00275                  const Double& sinDPA, const Double& cosDPA,
00276                  const Vector<Int>& Origin, const Vector<Int>& size);
00277 
00278     Complex* getConvFunc_p(Vector<Int>& cfShape,
00279                            CFBuffer& cfb,
00280                            Double& wVal, Int& fndx, 
00281                            Int& wndx,
00282                            PolMapType& mNdx, PolMapType& conjMNdx,
00283                            Int& ipol, uInt& mRow);
00284     void cachePhaseGrad_p(const Vector<Double>& pointingOffset,
00285                           const Vector<Int>&cfShape,
00286                           const Vector<Int>& convOrigin,
00287                           const Double& cfRefFreq,
00288                           const Double& imRefFreq,
00289                           const Int& spwID=0, const Int& fieldId=0);
00290   };
00291 }; //# NAMESPACE CASA - END
00292 
00293 #endif //