Line data Source code
1 : // -*- C++ -*-
2 : //# AWVisResampler.h: Definition of the AWVisResampler class
3 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 :
29 : #ifndef SYNTHESIS_AWVISRESAMPLER_H
30 : #define SYNTHESIS_AWVISRESAMPLER_H
31 :
32 : #include <synthesis/TransformMachines/CFStore.h>
33 : #include <synthesis/TransformMachines/VBStore.h>
34 : #include <synthesis/TransformMachines/VisibilityResampler.h>
35 : #include <msvis/MSVis/VisBuffer.h>
36 : #include <casacore/casa/Arrays/Array.h>
37 : #include <casacore/casa/Arrays/Vector.h>
38 :
39 : #include <casacore/casa/Logging/LogIO.h>
40 : #include <casacore/casa/Logging/LogSink.h>
41 : #include <casacore/casa/Logging/LogMessage.h>
42 :
43 : namespace casa { //# NAMESPACE CASA - BEGIN
44 : using namespace casacore;
45 : class AWVisResampler: public VisibilityResampler
46 : {
47 : public:
48 0 : AWVisResampler(): VisibilityResampler(),
49 : cached_phaseGrad_p(),
50 0 : cached_PointingOffset_p()
51 0 : {cached_PointingOffset_p.resize(2);cached_PointingOffset_p=-1000.0;runTimeG_p=runTimeDG_p=0.0;};
52 : // AWVisResampler(const CFStore& cfs): VisibilityResampler(cfs) {}
53 0 : virtual ~AWVisResampler() {};
54 :
55 0 : virtual VisibilityResamplerBase* clone()
56 0 : {return new AWVisResampler(*this);}
57 :
58 : // AWVisResampler(const AWVisResampler& other): VisibilityResampler(other),cfMap_p(), conjCFMap_p()
59 : // {copy(other);}
60 :
61 0 : virtual void copyMaps(const AWVisResampler& other)
62 0 : {setCFMaps(other.cfMap_p, other.conjCFMap_p);}
63 0 : virtual void copy(const VisibilityResamplerBase& other)
64 : {
65 0 : VisibilityResampler::copy(other);
66 : // const casacore::Vector<casacore::Int> cfmap=other.getCFMap();
67 : // const casacore::Vector<casacore::Int> conjcfmap = other.getConjCFMap();
68 :
69 : // setCFMaps(cfmap,conjcfmap);
70 0 : }
71 :
72 0 : virtual void copy(const AWVisResampler& other)
73 : {
74 0 : VisibilityResampler::copy(other);
75 0 : SynthesisUtils::SETVEC(cached_phaseGrad_p, other.cached_phaseGrad_p);
76 0 : SynthesisUtils::SETVEC(cached_PointingOffset_p, other.cached_PointingOffset_p);
77 0 : }
78 :
79 : AWVisResampler& operator=(const AWVisResampler& other)
80 : {
81 : copy(other);
82 : SynthesisUtils::SETVEC(cached_phaseGrad_p, other.cached_phaseGrad_p);
83 : SynthesisUtils::SETVEC(cached_PointingOffset_p, other.cached_PointingOffset_p);
84 : return *this;
85 : }
86 :
87 0 : virtual void setCFMaps(const casacore::Vector<casacore::Int>& cfMap, const casacore::Vector<casacore::Int>& conjCFMap)
88 0 : {SETVEC(cfMap_p,cfMap);SETVEC(conjCFMap_p,conjCFMap);}
89 :
90 : // virtual void setConvFunc(const CFStore& cfs) {convFuncStore_p = cfs;};
91 : //
92 : //------------------------------------------------------------------------------
93 : //
94 : // Re-sample the griddedData on the VisBuffer (a.k.a gridding).
95 : //
96 : // In this class, these just call the private templated version.
97 : // The first variant grids onto a double precision grid while the
98 : // second one does it on a single precision grid.
99 : //
100 : // Note that the following calls allow using any CFStore object
101 : // for gridding while de-gridding uses the internal
102 : // convFuncStore_p object.
103 : // virtual void DataToGrid(casacore::Array<casacore::DComplex>& griddedData, VBStore& vbs, casacore::Matrix<casacore::Double>& sumwt,
104 : // const casacore::Bool& dopsf, CFStore& cfs)
105 : // {DataToGridImpl_p(griddedData, vbs, sumwt,dopsf,cfs);}
106 :
107 : // virtual void DataToGrid(casacore::Array<casacore::Complex>& griddedData, VBStore& vbs, casacore::Matrix<casacore::Double>& sumwt,
108 : // const casacore::Bool& dopsf, CFStore& cfs)
109 : // {DataToGridImpl_p(griddedData, vbs, sumwt,dopsf,cfs);}
110 : //
111 : // Simulating defaulting CFStore arguemnt in the above calls to convFuncStore_p
112 : //
113 :
114 : //***TEMP REMOVAL OF casacore::DComplex gridder*****
115 :
116 0 : virtual void DataToGrid(casacore::Array<casacore::DComplex>& griddedData, VBStore& vbs, casacore::Matrix<casacore::Double>& sumwt,
117 : const casacore::Bool& dopsf,casacore::Bool useConjFreqCF=false)
118 0 : {DataToGridImpl_p(griddedData, vbs, sumwt,dopsf,useConjFreqCF);}
119 :
120 0 : virtual void DataToGrid(casacore::Array<casacore::Complex>& griddedData, VBStore& vbs, casacore::Matrix<casacore::Double>& sumwt,
121 : const casacore::Bool& dopsf,casacore::Bool useConjFreqCF=false)
122 0 : {DataToGridImpl_p(griddedData, vbs, sumwt,dopsf,useConjFreqCF);}
123 :
124 : //
125 : //------------------------------------------------------------------------------
126 : //
127 : // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding)
128 : //
129 : virtual void GridToData(VBStore& vbs,const casacore::Array<casacore::Complex>& griddedData);
130 : // virtual void GridToData(VBStore& vbs, casacore::Array<casacore::Complex>& griddedData);
131 : protected:
132 0 : virtual casacore::Complex getConvFuncVal(const casacore::Cube<casacore::Double>& convFunc, const casacore::Matrix<casacore::Double>& uvw,
133 : const casacore::Int& irow, const casacore::Vector<casacore::Int>& pixel)
134 : {
135 0 : (void)uvw; (void)irow;return convFunc(pixel[0],pixel[1],pixel[2]);
136 : }
137 : casacore::Complex getCFArea(casacore::Complex* __restrict__& convFuncV, casacore::Double& wVal,
138 : casacore::Vector<casacore::Int>& scaledSupport, casacore::Vector<casacore::Float>& scaledSampling,
139 : casacore::Vector<casacore::Double>& off,
140 : casacore::Vector<casacore::Int>& convOrigin, casacore::Vector<casacore::Int>& cfShape,
141 : casacore::Double& sinDPA, casacore::Double& cosDPA);
142 :
143 : template <class T>
144 : casacore::Complex accumulateOnGrid(casacore::Array<T>& grid, casacore::Complex* __restrict__& convFuncV,
145 : casacore::Complex& nvalue,
146 : casacore::Double& wVal, casacore::Vector<casacore::Int>& scaledSupport,
147 : casacore::Vector<casacore::Float>& scaledSampling, casacore::Vector<casacore::Double>& off,
148 : casacore::Vector<casacore::Int>& convOrigin, casacore::Vector<casacore::Int>& /*cfShape*/,
149 : casacore::Vector<casacore::Int>& loc, casacore::Vector<casacore::Int>& igrdpos,
150 : casacore::Double& /*sinDPA*/, casacore::Double& /*cosDPA*/,
151 : casacore::Bool& finitePointingOffset, casacore::Bool dopsf);
152 : template <class T>
153 : void XInnerLoop(const casacore::Int *scaleSupport, const casacore::Float* scaledSampling,
154 : const casacore::Double* off,
155 : const casacore::Int* loc, casacore::Complex& cfArea,
156 : const casacore::Int * __restrict__ iGrdPosPtr,
157 : casacore::Complex *__restrict__& convFuncV,
158 : const casacore::Int* convOrigin,
159 : casacore::Complex& nvalue,
160 : casacore::Double& wVal,
161 : casacore::Bool& /*finitePointingOffset*/,
162 : casacore::Bool& /*doPSFOnly*/,
163 : T* __restrict__ gridStore,
164 : casacore::Int* iloc,
165 : casacore::Complex& norm,
166 : casacore::Int* igrdpos);
167 :
168 : template <class T>
169 : void accumulateFromGrid(T& nvalue, const T* __restrict__& grid,
170 : casacore::Vector<casacore::Int>& iGrdPos,
171 : casacore::Complex* __restrict__& convFuncV,
172 : casacore::Double& wVal, casacore::Vector<casacore::Int>& scaledSupport,
173 : casacore::Vector<casacore::Float>& scaledSampling, casacore::Vector<casacore::Double>& off,
174 : casacore::Vector<casacore::Int>& convOrigin, casacore::Vector<casacore::Int>& cfShape,
175 : casacore::Vector<casacore::Int>& loc,
176 : casacore::Complex& phasor,
177 : casacore::Double& sinDPA, casacore::Double& cosDPA,
178 : casacore::Bool& finitePointingOffset,
179 : casacore::Matrix<casacore::Complex>& cached_phaseGrad_p);
180 :
181 : //
182 : //------------------------------------------------------------------------------
183 : //----------------------------Private parts-------------------------------------
184 : //------------------------------------------------------------------------------
185 : //
186 : private:
187 : // casacore::Vector<casacore::Double> uvwScale_p, offset_p, dphase_p;
188 : // casacore::Vector<casacore::Int> chanMap_p, polMap_p;
189 : // CFStore convFuncStore_p;
190 : // // casacore::Int inc0_p, inc1_p, inc2_p, inc3_p;
191 : // casacore::Vector<casacore::Int> inc_p;
192 : // casacore::Vector<casacore::Int> cfMap_p, conjCFMap_p;
193 : casacore::Vector<casacore::Int> gridInc_p, cfInc_p;
194 : casacore::Matrix<casacore::Complex> cached_phaseGrad_p;
195 : casacore::Vector<casacore::Double> cached_PointingOffset_p;
196 : //
197 : // Re-sample the griddedData on the VisBuffer (a.k.a de-gridding).
198 : //
199 : template <class T>
200 : void DataToGridImpl_p(casacore::Array<T>& griddedData, VBStore& vb,
201 : casacore::Matrix<casacore::Double>& sumwt,const casacore::Bool& dopsf,
202 : casacore::Bool /*useConjFreqCF*/);
203 :
204 : void sgrid(casacore::Vector<casacore::Double>& pos, casacore::Vector<casacore::Int>& loc, casacore::Vector<casacore::Double>& off,
205 : casacore::Complex& phasor, const casacore::Int& irow, const casacore::Matrix<casacore::Double>& uvw,
206 : const casacore::Double& dphase, const casacore::Double& freq,
207 : const casacore::Vector<casacore::Double>& scale, const casacore::Vector<casacore::Double>& offset,
208 : const casacore::Vector<casacore::Float>& sampling);
209 :
210 0 : inline casacore::Bool onGrid (const casacore::Int& nx, const casacore::Int& ny, const casacore::Int& nw,
211 : const casacore::Vector<casacore::Int>& loc,
212 : const casacore::Vector<casacore::Int>& support)
213 : {
214 0 : return (((loc(0)-support[0]) >= 0 ) && ((loc(0)+support[0]) < nx) &&
215 0 : ((loc(1)-support[1]) >= 0 ) && ((loc(1)+support[1]) < ny) &&
216 0 : (loc(2) >= 0) && (loc(2) <= nw));
217 : };
218 :
219 : // casacore::Array assignment operator in CASACore requires lhs.nelements()
220 : // == 0 or lhs.nelements()=rhs.nelements()
221 : template <class T>
222 0 : inline void SETVEC(casacore::Vector<T>& lhs, const casacore::Vector<T>& rhs)
223 0 : {lhs.resize(rhs.shape()); lhs = rhs;};
224 :
225 :
226 : //
227 : // Internal methods to address a 4D array. These should ulimately
228 : // moved to a Array4D class in CASACore
229 : //
230 :
231 : // This is called less frequently. Currently once per VisBuffer
232 : // inline void cacheAxisIncrements(const casacore::Vector<casacore::Int>& n, casacore::Vector<casacore::Int>& inc)
233 : // {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];}
234 :
235 :
236 : // The following method is also called from the inner loop, but
237 : // does not use CASA casacore::Vector (which are not thread safe, I (SB) am
238 : // told).
239 0 : inline casacore::Complex getFrom4DArray(const casacore::Complex *__restrict__& store,
240 : const casacore::Int* iPos, const casacore::Int* inc)
241 : {
242 0 : return *(store+(iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]));
243 : // return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];
244 : };
245 :
246 : // The following two methods are called in the innermost loop.
247 0 : inline casacore::Complex getFrom4DArray(const casacore::Complex *__restrict__& store,
248 : const casacore::Vector<casacore::Int>& iPos, const casacore::Vector<casacore::Int>& inc)
249 : {
250 0 : return *(store+(iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]));
251 : // return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];
252 : };
253 : inline casacore::DComplex getFrom4DArray(const casacore::DComplex *__restrict__& store,
254 : const casacore::Vector<casacore::Int>& iPos, const casacore::Vector<casacore::Int>& inc)
255 : {
256 : return *(store+(iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]));
257 : // return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];
258 : };
259 :
260 : template <class T>
261 0 : void addTo4DArray(T *__restrict__& store,
262 : const casacore::Int *__restrict__& iPos, const casacore::Vector<casacore::Int>& inc,
263 : casacore::Complex& nvalue, casacore::Complex& wt) __restrict__
264 : {
265 : // T *tmp=store+(iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]);
266 : // *tmp += nvalue*wt;
267 0 : store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]] += (nvalue*wt);
268 0 : }
269 :
270 : //
271 : // This rotates the convolution function by rotating the
272 : // co-ordinate system. For the accuracies already required for
273 : // EVLA and ALMA, this is not useful. Leaving it hear for now....
274 : //
275 : casacore::Bool reindex(const casacore::Vector<casacore::Int>& in, casacore::Vector<casacore::Int>& out,
276 : const casacore::Double& sinDPA, const casacore::Double& cosDPA,
277 : const casacore::Vector<casacore::Int>& Origin, const casacore::Vector<casacore::Int>& size);
278 :
279 : casacore::Complex* getConvFunc_p(const double& vbPA,
280 : casacore::Vector<casacore::Int>& cfShape,
281 : casacore::Vector<int>& support,
282 : int& muellerElement,
283 : CFBuffer& cfb,
284 : casacore::Double& wVal, casacore::Int& fndx,
285 : casacore::Int& wndx,
286 : PolMapType& mNdx, PolMapType& conjMNdx,
287 : casacore::Int& ipol, casacore::uInt& mRow);
288 : void cachePhaseGrad_p(const casacore::Vector<casacore::Double>& pointingOffset,
289 : const casacore::Vector<casacore::Int>&cfShape,
290 : const casacore::Vector<casacore::Int>& convOrigin,
291 : const casacore::Double& cfRefFreq,
292 : const casacore::Double& imRefFreq,
293 : const casacore::Int& spwID=0, const casacore::Int& fieldId=0);
294 : };
295 : }; //# NAMESPACE CASA - END
296 :
297 : #endif //
|
