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