Line data Source code
1 : // -*- C++ -*-
2 : //# AWVisResampler.cc: Implementation 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: casa-feedback@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 : #include <synthesis/TransformMachines/SynthesisError.h>
30 : #include <synthesis/TransformMachines/AWVisResampler.h>
31 : #include <synthesis/TransformMachines/Utils.h>
32 : #include <synthesis/TransformMachines/SynthesisMath.h>
33 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
34 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
35 : #include <casacore/casa/OS/Timer.h>
36 : #include <fstream>
37 : #include <iostream>
38 : #include <typeinfo>
39 : #include <iomanip>
40 : #include <synthesis/TransformMachines/FortranizedLoops.h>
41 : #ifdef _OPENMP
42 : #include <omp.h>
43 : #endif
44 : extern "C" {
45 : void clLoopsToGrid();
46 : };
47 : //#include <casa/BasicMath/Functors.h>
48 : using namespace casacore;
49 : namespace casa{
50 :
51 : //
52 : //-----------------------------------------------------------------------------------
53 : // Re-sample the griddedData on the VisBuffer (a.k.a gridding)
54 : //
55 : // Template instantiations for re-sampling onto a double precision
56 : // or single precision grid.
57 : //
58 : template
59 : void AWVisResampler::addTo4DArray(DComplex* __restrict__ & store,
60 : const Int* __restrict__ & iPos,
61 : const Vector<Int>& inc,
62 : Complex& nvalue, Complex& wt) __restrict__ ;
63 : template
64 : void AWVisResampler::addTo4DArray(Complex* __restrict__ & store,
65 : const Int* __restrict__ & iPos,
66 : const Vector<Int>& inc,
67 : Complex& nvalue, Complex& wt) __restrict__;
68 : template
69 : void AWVisResampler::DataToGridImpl_p(Array<DComplex>& grid, VBStore& vbs,
70 : Matrix<Double>& sumwt,const Bool& dopsf,
71 : Bool useConjFreqCF); // __restrict__;
72 : template
73 : void AWVisResampler::DataToGridImpl_p(Array<Complex>& grid, VBStore& vbs,
74 : Matrix<Double>& sumwt,const Bool& dopsf,
75 : Bool useConjFreqCF); // __restrict__;
76 :
77 : template
78 : Complex AWVisResampler::accumulateOnGrid(Array<DComplex>& grid, Complex* __restrict__& convFuncV,
79 : Complex& nvalue, Double& wVal,
80 : Vector<Int>& scaledSupport, Vector<Float>& scaledSampling,
81 : Vector<Double>& off, Vector<Int>& convOrigin,
82 : Vector<Int>& /*cfShape*/, Vector<Int>& loc, Vector<Int>& igrdpos,
83 : Double& /*sinDPA*/, Double& /*cosDPA*/,
84 : Bool& pointingOffset, Bool dopsf);
85 : template
86 : Complex AWVisResampler::accumulateOnGrid(Array<Complex>& grid, Complex* __restrict__& convFuncV,
87 : Complex& nvalue, Double& wVal,
88 : Vector<Int>& scaledSupport, Vector<Float>& scaledSampling,
89 : Vector<Double>& off, Vector<Int>& convOrigin,
90 : Vector<Int>& /*cfShape*/, Vector<Int>& loc, Vector<Int>& igrdpos,
91 : Double& /*sinDPA*/, Double& /*cosDPA*/,
92 : Bool& pointingOffset, Bool dopsf);
93 : // template
94 : // void AWVisResampler::accumulateFromGrid(Complex& nvalue, const DComplex* __restrict__& grid,
95 : // Vector<Int>& iGrdPos,
96 : // Complex* __restrict__& convFuncV,
97 : // Double& wVal, Vector<Int>& scaledSupport,
98 : // Vector<Float>& scaledSampling, Vector<Double>& off,
99 : // Vector<Int>& convOrigin, Vector<Int>& cfShape,
100 : // Vector<Int>& loc,
101 : // Complex& phasor,
102 : // Double& sinDPA, Double& cosDPA,
103 : // Bool& finitePointingOffset,
104 : // Matrix<Complex>& cached_phaseGrad_p,
105 : // Bool dopsf);
106 : template
107 : void AWVisResampler::accumulateFromGrid(Complex& nvalue, const Complex* __restrict__& grid,
108 : Vector<Int>& iGrdPos,
109 : Complex* __restrict__& convFuncV,
110 : Double& wVal, Vector<Int>& scaledSupport,
111 : Vector<Float>& scaledSampling, Vector<Double>& off,
112 : Vector<Int>& convOrigin, Vector<Int>& cfShape,
113 : Vector<Int>& loc,
114 : Complex& phasor,
115 : Double& sinDPA, Double& cosDPA,
116 : Bool& finitePointingOffset,
117 : Matrix<Complex>& cached_phaseGrad_p);
118 :
119 : template
120 : void AWVisResampler::XInnerLoop(const Int *scaledSupport, const Float* scaledSampling,
121 : const Double* off,
122 : const Int* loc, Complex& cfArea,
123 : const Int * __restrict__ iGrdPosPtr,
124 : Complex *__restrict__& convFuncV,
125 : const Int* convOrigin,
126 : Complex& nvalue,
127 : Double& wVal,
128 : Bool& /*finitePointingOffset*/,
129 : Bool& /*doPSFOnly*/,
130 : Complex* __restrict__ gridStore,
131 : Int* iloc,
132 : Complex& norm,
133 : Int* igrdpos);
134 : template
135 : void AWVisResampler::XInnerLoop(const Int *scaledSupport, const Float* scaledSampling,
136 : const Double* off,
137 : const Int* loc, Complex& cfArea,
138 : const Int * __restrict__ iGrdPosPtr,
139 : Complex *__restrict__& convFuncV,
140 : const Int* convOrigin,
141 : Complex& nvalue,
142 : Double& wVal,
143 : Bool& /*finitePointingOffset*/,
144 : Bool& /*doPSFOnly*/,
145 : DComplex* __restrict__ gridStore,
146 : Int* iloc,
147 : Complex& norm,
148 : Int* igrdpos);
149 :
150 0 : Complex* AWVisResampler::getConvFunc_p(const double& vbPA, Vector<Int>& cfShape,
151 : Vector<int>& support,
152 : int& muellerElement,
153 : CFBuffer& cfb,
154 : Double& wVal, Int& fndx, Int& wndx,
155 : PolMapType& mNdx, PolMapType& conjMNdx,
156 : Int& ipol, uInt& mRow)
157 : {
158 : Bool Dummy;
159 : Array<Complex> *convFuncV;
160 : CFCell *cfcell;
161 : //
162 : // Since we conjugate the CF depending on the sign of the w-value,
163 : // pick the appropriate element of the Mueller Matrix (see note on
164 : // this for details). Without reading the note/understanding,
165 : // fiddle with this logic at your own risk (can easily lead to a
166 : // lot of grief. --Sanjay).
167 : //
168 0 : timer_p.mark();
169 :
170 0 : if (wVal > 0.0)
171 : {
172 0 : cfcell=&(*(cfb.getCFCellPtr(fndx,wndx,mNdx[ipol][mRow])));
173 0 : CFCell& cfO=cfb(fndx,wndx,mNdx[ipol][mRow]);
174 0 : convFuncV = &(*cfO.getStorage());
175 0 : support(0)=support(1)=cfO.xSupport_p;
176 : // convFuncV=&(*(cfb.getCFCellPtr(fndx,wndx,mNdx[ipol][mRow])->getStorage()));//->getStorage(Dummy);
177 : // if (mNdx[ipol][mRow] != ipol)
178 : // cerr << "Indexes+: " << fndx << " " << wndx << " " << mNdx[ipol][mRow] << " " << ipol << " " << mRow << endl;
179 : }
180 : else
181 : {
182 0 : cfcell=&(*(cfb.getCFCellPtr(fndx,wndx,conjMNdx[ipol][mRow])));
183 0 : CFCell& cfO=cfb(fndx,wndx,conjMNdx[ipol][mRow]);
184 0 : convFuncV = &(*cfO.getStorage());
185 0 : support(0)=support(1)=cfO.xSupport_p;
186 : // convFuncV=&(*(cfb.getCFCellPtr(fndx,wndx,conjMNdx[ipol][mRow])->getStorage()));//->getStorage(Dummy);
187 : // if (conjMNdx[ipol][mRow] != ipol)
188 : // cerr << "Indexes-: " << fndx << " " << wndx << " " << conjMNdx[ipol][mRow] << " " << ipol << " " << mRow << endl;
189 : }
190 :
191 : // cerr << getpid() << " " << cfb.getCFCacheDir() << " " << cfcell->fileName_p << " " << cfcell->freqValue_p << endl;
192 :
193 :
194 : // Get the pointer to the CFCell storage (a single CF)
195 : // if ((convFuncV = &(*cfcell->getStorage())) == NULL)
196 0 : if (convFuncV == NULL)
197 0 : throw(SynthesisFTMachineError("cfcell->getStorage() == null"));
198 :
199 : // Load the CF if it not already loaded. If a new CF is loaded,
200 : // check if it needs to be rotated.
201 0 : if (convFuncV->shape().product() == 0)
202 : {
203 0 : Array<Complex> tt=SynthesisUtils::getCFPixels(cfb.getCFCacheDir(), cfcell->fileName_p);
204 0 : cfcell->setStorage(tt);
205 :
206 : //cerr << (cfcell->isRotationallySymmetric_p?"o":"+");
207 :
208 : // No rotation necessary if the CF is rotationally symmetric
209 0 : if (!(cfcell->isRotationallySymmetric_p))
210 : {
211 0 : CFCell *baseCFC=NULL;
212 : // Rotate only if the difference between CF PA and VB PA
213 : // is greater than paTolerance.
214 0 : SynthesisUtils::rotate2(vbPA, *baseCFC, *cfcell, paTolerance_p);
215 : }
216 0 : convFuncV = &(*cfcell->getStorage());
217 0 : }
218 :
219 : //cfShape.reference(cfcell->cfShape_p);
220 0 : cfShape.assign(convFuncV->shape().asVector());
221 :
222 : // Always extract the Mueller element value from mNdx. mNdx
223 : // carries the direct mapping between Mueller Matrix and
224 : // Visibility vector.
225 0 : muellerElement=cfb.getCFCellPtr(fndx,wndx,mNdx[ipol][mRow])->muellerElement_p;
226 :
227 0 : runTimeG1_p += timer_p.real();
228 :
229 :
230 0 : return convFuncV->getStorage(Dummy);
231 : };
232 :
233 : template <class T>
234 0 : void AWVisResampler::XInnerLoop(const Int *scaledSupport, const Float* scaledSampling,
235 : const Double* off,
236 : const Int* loc, Complex& cfArea,
237 : const Int * __restrict__ iGrdPosPtr,
238 : Complex *__restrict__& convFuncV,
239 : const Int* convOrigin,
240 : Complex& nvalue,
241 : Double& wVal,
242 : Bool& /*finitePointingOffset*/,
243 : Bool& /*doPSFOnly*/,
244 : T* __restrict__ gridStore,
245 : Int* iloc,
246 : Complex& norm,
247 : Int* igrdpos)
248 : {
249 0 : Complex wt;
250 0 : const Int *tt=iloc;
251 : Bool Dummy;
252 0 : const Int *cfInc_ptr=cfInc_p.getStorage(Dummy);
253 0 : for(Int ix=-scaledSupport[0]; ix <= scaledSupport[0]; ix++)
254 : {
255 0 : iloc[0]=(Int)((scaledSampling[0]*ix+off[0])-1)+convOrigin[0];
256 0 : igrdpos[0]=loc[0]+ix;
257 :
258 : {
259 0 : wt = getFrom4DArray((const Complex * __restrict__ &)convFuncV,
260 0 : tt,cfInc_ptr)/cfArea;
261 0 : if (wVal > 0.0) {wt = conj(wt);}
262 0 : norm += (wt);
263 : // if (finitePointingOffset && !doPSFOnly)
264 : // wt *= cached_phaseGrad_p(iloc[0]+phaseGradOrigin_l[0],
265 : // iloc[1]+phaseGradOrigin_l[1]);
266 :
267 : // The following uses raw index on the 4D grid
268 0 : addTo4DArray(gridStore,iGrdPosPtr,gridInc_p, nvalue,wt);
269 : }
270 : }
271 0 : }
272 :
273 : template <class T>
274 0 : Complex AWVisResampler::accumulateOnGrid(Array<T>& grid,Complex* __restrict__& convFuncV,
275 : Complex& nvalue,Double& wVal,
276 : Vector<Int>& scaledSupport, Vector<Float>& scaledSampling,
277 : Vector<Double>& off, Vector<Int>& convOrigin,
278 : Vector<Int>& /*cfShape*/, Vector<Int>& loc, Vector<Int>& igrdpos,
279 : Double& /*sinDPA*/, Double& /*cosDPA*/,
280 : Bool& finitePointingOffset,
281 : Bool doPSFOnly)
282 : {
283 0 : Vector<Int> iloc(4,0), tiloc(4);
284 : Bool Dummy;
285 : // wt no longer appears to be used
286 : // Complex wt
287 0 : Complex cfArea=1.0;
288 0 : Complex norm=0.0;
289 0 : const Int * __restrict__ iGrdPosPtr = igrdpos.getStorage(Dummy);
290 0 : T* __restrict__ gridStore = grid.getStorage(Dummy);
291 : // Nth no longer appears to be used
292 : // Int Nth = 1;
293 : //#ifdef _OPENMP
294 : // Nth=max(omp_get_max_threads()-2,1);
295 : //#endif
296 :
297 0 : const Int* scaledSupport_ptr=scaledSupport.getStorage(Dummy);
298 0 : const Float *scaledSampling_ptr=scaledSampling.getStorage(Dummy);
299 0 : const Double *off_ptr=off.getStorage(Dummy);
300 0 : const Int *loc_ptr = loc.getStorage(Dummy);
301 0 : const Int* convOrigin_ptr=convOrigin.getStorage(Dummy);
302 0 : Int *iloc_ptr=iloc.getStorage(Dummy);
303 0 : Int *igrdpos_ptr=igrdpos.getStorage(Dummy);
304 :
305 0 : Bool finitePointingOffset_l=finitePointingOffset;
306 0 : Bool doPSFOnly_l=doPSFOnly;
307 0 : Double wVal_l=wVal;
308 0 : Complex nvalue_l=nvalue;
309 0 : Complex *convFuncV_l=convFuncV;
310 : //
311 : // !!!!! Compute cfArea for high precision work
312 : //
313 :
314 : // cfArea = getCFArea(convFuncV, wVal, scaledSupport, scaledSampling,
315 : // off, convOrigin, cfShape, sinDPA,cosDPA);
316 :
317 :
318 : // cerr << "Cfarea = " << cfArea << endl;
319 0 : IPosition phaseGradOrigin_l;
320 0 : phaseGradOrigin_l = cached_phaseGrad_p.shape()/2;
321 :
322 : // #pragma omp parallel default(none) \
323 : // shared(gridStore) \
324 : // firstprivate(scaledSupport_ptr,scaledSampling_ptr,off_ptr,loc_ptr,cfArea,iGrdPosPtr, \
325 : // convFuncV_l, convOrigin_ptr, nvalue_l, wVal_l, finitePointingOffset_l, doPSFOnly_l, \
326 : // iloc_ptr, norm,igrdpos_ptr) num_threads(Nth)
327 : {
328 : // #pragma omp for
329 0 : for(Int iy=-scaledSupport[1]; iy <= scaledSupport[1]; iy++)
330 : {
331 0 : iloc_ptr[1]=(Int)((scaledSampling[1]*iy+off[1])-1)+convOrigin[1];
332 0 : igrdpos[1]=loc[1]+iy;
333 0 : XInnerLoop(scaledSupport_ptr, scaledSampling_ptr,
334 : off_ptr,
335 : loc_ptr, cfArea,
336 : iGrdPosPtr,
337 : convFuncV_l,
338 : convOrigin_ptr,
339 : nvalue_l,
340 : wVal_l,
341 : finitePointingOffset_l,
342 : doPSFOnly_l,
343 : gridStore,
344 : iloc_ptr,
345 : norm,
346 : igrdpos_ptr);
347 :
348 : // for(Int ix=-scaledSupport[0]; ix <= scaledSupport[0]; ix++)
349 : // {
350 : // iloc[0]=(Int)((scaledSampling[0]*ix+off[0])-1)+convOrigin[0];
351 : // igrdpos[0]=loc[0]+ix;
352 : // //
353 : // // reindex() is for historical reasons and does three
354 : // // operations: (1) rotate the co-ord. sys. using
355 : // // sin/cosDPA, (2) add convOrigin to iloc and return the
356 : // // result in tiloc and add convOrigin to tiloc, and (3)
357 : // // return true if tiloc lines with in the cfShape.
358 : // //
359 : // // if (reindex(iloc,tiloc,sinDPA, cosDPA, convOrigin, cfShape))
360 : // {
361 : // wt = getFrom4DArray((const Complex * __restrict__ &)convFuncV,
362 : // iloc,cfInc_p)/cfArea;
363 : // if (wVal > 0.0) {wt = conj(wt);}
364 : // norm += (wt);
365 : // if (finitePointingOffset && !doPSFOnly)
366 : // wt *= cached_phaseGrad_p(iloc[0]+phaseGradOrigin_l[0],
367 : // iloc[1]+phaseGradOrigin_l[1]);
368 : // // The following uses raw index on the 4D grid
369 : // addTo4DArray(gridStore,iGrdPosPtr,gridInc_p, nvalue,wt);
370 : // }
371 : // }
372 : }
373 : }
374 0 : return norm;
375 0 : }
376 : // Moved the accumulateFromGrid() method to file to play with
377 : // multi-threading it to not clutter this file. Both versions
378 : // (threaded and non-threaded) are in this file.
379 : #include "accumulateFromGrid.inc"
380 : //
381 : //-----------------------------------------------------------------------------------
382 : //
383 0 : void AWVisResampler::cachePhaseGrad_p(const Vector<Double>& pointingOffset,
384 : const Vector<Int>&cfShape,
385 : const Vector<Int>& convOrigin,
386 : const Double& /*cfRefFreq*/,
387 : const Double& /*imRefFreq*/,
388 : const Int& spwID, const Int& fieldId)
389 : {
390 0 : LogIO log_l(LogOrigin("AWVisResampler","cachePhaseGrad[R&D]"));
391 : //cout << "# " << cfRefFreq << " " << imRefFreq << endl;
392 0 : if (
393 0 : (sum(fabs(pointingOffset-cached_PointingOffset_p)) > 1e-6) ||
394 0 : (cached_phaseGrad_p.shape()[0] < cfShape[0]) ||
395 0 : (cached_phaseGrad_p.shape()[1] < cfShape[1])
396 : )
397 : {
398 : log_l << "Computing phase gradiant for pointing offset "
399 : << pointingOffset << cfShape << " " << cached_phaseGrad_p.shape()
400 : << "(SPW: " << spwID << " Field: " << fieldId << ")"
401 0 : << LogIO::DEBUGGING << LogIO::POST;
402 0 : Int nx=cfShape(0), ny=cfShape(1);
403 : Double grad;
404 0 : Complex phx,phy;
405 :
406 0 : cached_phaseGrad_p.resize(nx,ny);
407 0 : cached_PointingOffset_p = pointingOffset;
408 :
409 0 : for(Int ix=0;ix<nx;ix++)
410 : {
411 0 : grad = (ix-convOrigin[0])*pointingOffset[0];
412 : Double sx,cx;
413 0 : SINCOS(grad,sx,cx);
414 : // phx = Complex(cos(grad),sin(grad));
415 0 : phx = Complex(cx,sx);
416 0 : for(Int iy=0;iy<ny;iy++)
417 : {
418 0 : grad = (iy-convOrigin[1])*pointingOffset[1];
419 : Double sy,cy;
420 0 : SINCOS(grad,sy,cy);
421 : // phy = Complex(cos(grad),sin(grad));
422 0 : phy = Complex(cy,sy);
423 0 : cached_phaseGrad_p(ix,iy)=phx*phy;
424 : }
425 : }
426 : }
427 0 : }
428 : //
429 : //-----------------------------------------------------------------------------------
430 : //
431 : // AWVisResampler& AWVisResampler::operator=(const AWVisResampler& other)
432 : // {
433 : // SETVEC(uvwScale_p, other.uvwScale_p);
434 : // SETVEC(offset_p, other.offset_p);
435 : // SETVEC(dphase_p, other.dphase_p);
436 : // SETVEC(chanMap_p, other.chanMap_p);
437 : // SETVEC(polMap_p, other.polMap_p);
438 :
439 : // convFuncStore_p = other.convFuncStore_p;
440 :
441 : // return *this;
442 : // }
443 :
444 : // CFB::initPolMaps(polMap,conjPolMap) sets the internal maps of CFB
445 : //
446 : // AWP::findCF() --> cfCache->initPolMap(...,...) --> cfb->initPolMaps(...,...)
447 : //
448 : // AWVR extracts polMap from CFB.
449 : // Same CF is extracted for gridding and de-gridding. CF* used in the gridding loops.
450 : // getConvFunc_p() ensures the jugglery needed for AW CFs
451 :
452 : //
453 : //-----------------------------------------------------------------------------------
454 : // Template implementation for DataToGrid
455 : //
456 : template <class T>
457 0 : void AWVisResampler::DataToGridImpl_p(Array<T>& grid, VBStore& vbs,
458 : Matrix<Double>& sumwt,const Bool& dopsf,
459 : Bool /*useConjFreqCF*/)
460 : {
461 0 : LogIO log_l(LogOrigin("AWVisResampler[R&D]","DataToGridImpl_p"));
462 : Int nDataChan, nDataPol, nGridPol, nGridChan, nx, ny, nw;//, nCFFreq;
463 : Int targetIMChan, targetIMPol, rbeg, rend;//, PolnPlane, ConjPlane;
464 : Int startChan, endChan;
465 :
466 0 : Vector<Float> sampling(2),scaledSampling(2);
467 0 : Vector<Int> support(2),loc(3), iloc(4),tiloc(4),scaledSupport(2);
468 0 : Vector<Double> pos(3), off(3);
469 0 : Vector<Int> igrdpos(4);
470 :
471 0 : Complex phasor, nvalue, wt;
472 0 : Complex norm;
473 0 : Vector<Int> cfShape;
474 0 : cfShape=vbRow2CFBMap_p(0)->getStorage()(0,0,0)->getStorage()->shape().asVector();
475 :
476 0 : Vector<Int> convOrigin = (cfShape)/2;
477 0 : Double sinDPA=0.0, cosDPA=1.0, cfRefFreq;
478 : // Double cfScale=1.0;
479 :
480 0 : timer_p.mark();
481 0 : rbeg = 0; rend = vbs.nRow_p;
482 0 : rbeg = vbs.beginRow_p;
483 0 : rend = vbs.endRow_p;
484 :
485 0 : nx = grid.shape()[0]; ny = grid.shape()[1];
486 0 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
487 :
488 0 : nDataPol = vbs.flagCube_p.shape()[0];
489 0 : nDataChan = vbs.flagCube_p.shape()[1];
490 :
491 : Bool Dummy, gDummy,
492 0 : accumCFs=((vbs.uvw_p.nelements() == 0) && dopsf);
493 :
494 0 : T* __restrict__ gridStore = grid.getStorage(gDummy);
495 :
496 0 : Double *freq=vbs.freq_p.getStorage(Dummy);
497 :
498 0 : cacheAxisIncrements(grid.shape().asVector(), gridInc_p);
499 :
500 0 : Bool * __restrict__ flagCube_ptr=vbs.flagCube_p.getStorage(Dummy);
501 0 : Bool * __restrict__ rowFlag_ptr = vbs.rowFlag_p.getStorage(Dummy);;
502 0 : Float * __restrict__ imgWts_ptr = vbs.imagingWeight_p.getStorage(Dummy);
503 0 : Complex * __restrict__ visCube_ptr = vbs.visCube_p.getStorage(Dummy);
504 :
505 0 : Vector<Double> wVals, fVals; PolMapType mVals, mNdx, conjMVals, conjMNdx;
506 : Double fIncr, wIncr;
507 0 : CFBuffer& cfb = *vbRow2CFBMap_p(0);
508 : // CFBStruct cfbst;
509 : // cfb.getAsStruct(cfbst);
510 : // for(int ii=0;ii<vbs.cfBSt_p.shape[0];ii++)
511 : // for(int jj=0;jj<vbs.cfBSt_p.shape[1];jj++)
512 : // for(int kk=0;kk<vbs.cfBSt_p.shape[2];kk++)
513 : // {
514 : // CFCStruct cfcst=vbs.cfBSt_p.getCFB(ii,jj,kk);
515 : // cerr << "[" << ii << "," << jj << "," << kk << "]:"
516 : // << cfcst.sampling << " "
517 : // << cfcst.xSupport << " "
518 : // << cfcst.ySupport
519 : // << endl;
520 : // }
521 :
522 :
523 :
524 0 : cfb.getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
525 0 : Vector<Double> pointingOffset(cfb.getPointingOffset());
526 0 : runTimeG1_p += timer_p.real();
527 :
528 0 : nw = wVals.nelements();
529 : //nCFFreq = fVals.nelements()-1;
530 0 : iloc = 0;
531 :
532 : // timer.mark();
533 0 : IPosition shp=vbs.flagCube_p.shape();
534 0 : Cube<Bool> allPolNChanDone_l(shp(0),shp(1),1);
535 0 : if (accumCFs)
536 : {
537 0 : for (Int ipol=0;ipol<nDataPol;ipol++)
538 : {
539 0 : if (polMap_p(ipol) < 0)
540 : {
541 0 : for (Int ichan=0;ichan<nDataChan;ichan++)
542 : //for (Int irow=rbeg;irow<rend;irow++)
543 0 : allPolNChanDone_l(ipol,ichan,0)=true;
544 : }
545 : }
546 :
547 0 : startChan = vbs.startChan_p;
548 0 : endChan = vbs.endChan_p;
549 : }
550 : else
551 : {
552 0 : startChan = 0;
553 0 : endChan = nDataChan;
554 : }
555 :
556 0 : Bool finitePointingOffsets=(
557 0 : (fabs(pointingOffset(0))>0) ||
558 0 : (fabs(pointingOffset(1))>0)
559 : );
560 0 : Bool isGridSinglePrecision=(typeid(gridStore[0]) == typeid(wt));
561 0 : Int vbSpw = (vbs.vb_p)->spectralWindow();
562 :
563 : // cerr << "-------------------BEGIN VB-------------------------" << endl;
564 :
565 0 : for(Int irow=rbeg; irow< rend; irow++){
566 : // if ((vbs.uvw_p.nelements() == 0))
567 : //if (accumCFs) if (allTrue(allPolNChanDone_l)) break;
568 :
569 0 : if(!(*(rowFlag_ptr+irow)))
570 : {
571 0 : for(Int ichan=startChan; ichan< endChan; ichan++)
572 : {
573 0 : if (*(imgWts_ptr + ichan+irow*nDataChan)!=0.0)
574 : {
575 0 : targetIMChan=chanMap_p[ichan];
576 :
577 0 : if((targetIMChan>=0) && (targetIMChan<nGridChan))
578 : {
579 0 : timer_p.mark();
580 0 : Double dataWVal = vbs.vb_p->uvw()(irow)(2);
581 0 : Int wndx = cfb.nearestWNdx(abs(dataWVal)*freq[ichan]/C::c);
582 0 : Int cfFreqNdx = cfb.nearestFreqNdx(vbSpw,ichan,vbs.conjBeams_p);
583 0 : runTimeG3_p += timer_p.real();
584 :
585 : Float s;
586 : //
587 : //------------------------------------------------------------------------------
588 : //
589 : // Using the int-index version for Freq, W and Muellerelements
590 : // cfb.getParams(cfRefFreq, s, support(0), support(1),0,wndx,0);
591 : //
592 : //------------------------------------------------------------------------------
593 :
594 0 : timer_p.mark();
595 : //s=cfb(cfFreqNdx,wndx,0).sampling_p; //Sampling is the same for all pol. planes
596 0 : cfb.getParams(cfRefFreq, s, support(0), support(1),cfFreqNdx,wndx,0);
597 0 : sampling(0) = sampling(1) = SynthesisUtils::nint(s);
598 0 : runTimeG4_p += timer_p.real();
599 :
600 0 : timer_p.mark();
601 0 : sgrid(pos,loc,off, phasor, irow, vbs.uvw_p, dphase_p[irow], freq[ichan],
602 0 : uvwScale_p, offset_p, sampling);
603 0 : runTimeG5_p += timer_p.real();
604 :
605 : //if (onGrid(nx, ny, nw, loc, support))
606 : {
607 : // Loop over all image-plane polarization planes.
608 0 : for(Int ipol=0; ipol< nDataPol; ipol++)
609 : {
610 0 : if((!(*(flagCube_ptr + ipol + ichan*nDataPol + irow*nDataPol*nDataChan))))
611 : {
612 0 : targetIMPol=polMap_p(ipol);
613 0 : if ((targetIMPol>=0) && (targetIMPol<nGridPol))
614 : {
615 0 : igrdpos[2]=targetIMPol; igrdpos[3]=targetIMChan;
616 :
617 0 : if(accumCFs) allPolNChanDone_l(ipol,ichan,0)=true;
618 :
619 : // if(dopsf) nvalue=Complex(*(imgWts_ptr + ichan + irow*nDataChan));
620 : // else nvalue=Complex(*(imgWts_ptr+ichan+irow*nDataChan))*
621 : // (*(visCube_ptr+ipol+ichan*nDataPol+irow*nDataChan*nDataPol)*phasor);
622 :
623 0 : norm = 0.0;
624 : // Loop over all relevant elements of the Mueller matrix for the polarization
625 : // ipol.
626 0 : Vector<int> conjMRow = conjMNdx[ipol];
627 : //for (uInt mCols=0;mCols<conjMNdx[ipol].nelements(); mCols++)
628 :
629 : // ipol determines the targetIMPol. Each targetIMPol gets a row of CFs (mRow).
630 : // visVecElements is gridded using the convFuncV and added to the target grid.
631 :
632 0 : for (uInt mCols=0;mCols<conjMRow.nelements(); mCols++)
633 : {
634 0 : int visVecElement=mCols, muellerElement;
635 :
636 0 : Complex* convFuncV=NULL;
637 0 : timer_p.mark();
638 : try
639 : {
640 0 : convFuncV=getConvFunc_p(vbs.paQuant_p.getValue("rad"),
641 : cfShape, support,muellerElement,
642 : cfb, dataWVal, cfFreqNdx,
643 : wndx, mNdx, conjMNdx,
644 : ipol, mCols);
645 : }
646 0 : catch (SynthesisFTMachineError& x)
647 : {
648 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
649 : }
650 : // Extract the vis. vector element corresponding to the mCols column of the conjMRow row of the Mueller matrix.
651 :
652 0 : visVecElement=(int)(muellerElement%nDataPol);
653 : // If the vis. vector element is flagged, don't grid it.
654 0 : if(((*(flagCube_ptr + visVecElement + ichan*nDataPol + irow*nDataPol*nDataChan)))) break;
655 :
656 : //cerr << "G: " << mCols << "-->" << visVecElement << "-->" << ipol << " " << polMap_p[ipol] << endl;
657 :
658 0 : if(dopsf) nvalue=Complex(*(imgWts_ptr + ichan + irow*nDataChan));
659 0 : else nvalue=Complex(*(imgWts_ptr+ichan+irow*nDataChan))*
660 0 : (*(visCube_ptr+visVecElement+ichan*nDataPol+irow*nDataChan*nDataPol)*phasor);
661 :
662 0 : if (!onGrid(nx, ny, nw, loc, support)) break;
663 0 : runTimeG6_p += timer_p.real();
664 :
665 0 : convOrigin=cfShape/2;
666 0 : Bool psfOnly=((dopsf==true) && (accumCFs==false));
667 0 : if (finitePointingOffsets )
668 0 : cachePhaseGrad_p(pointingOffset, cfShape, convOrigin, cfRefFreq, vbs.imRefFreq(),
669 0 : ((const Int)(vbs.vb_p)->spectralWindow()),((const Int)((vbs.vb_p)->fieldId())));
670 :
671 0 : cacheAxisIncrements(cfShape, cfInc_p);
672 :
673 : // accumulateOnGrid() is a local C++ method with the inner loops. The include
674 : // file (FortanizedLoopsToGrid.cc) has the interface code to call the inner
675 : // loops re-written in FORTRAN (in synthesis/fortran/faccumulateOnGrid.f)
676 : #include <synthesis/TransformMachines/FortranizedLoopsToGrid.inc>
677 : }
678 0 : sumwt(targetIMPol,targetIMChan) += vbs.imagingWeight_p(ichan, irow)*abs(norm);
679 : // *(sumWt_ptr+apol+achan*nGridChan)+= *(imgWts_ptr+ichan+irow*nDataChan);
680 0 : }
681 : }
682 : } // End poln-loop
683 : } //End of on-grid loop
684 : }
685 : }
686 : } // End chan-loop
687 : }
688 : } // End row-loop
689 : //if (!dopsf) exit(0);
690 : // cerr << "-------------------END VB-------------------------" << endl;
691 :
692 0 : runTimeG_p = timer_p.real() + runTimeG1_p + runTimeG2_p + runTimeG3_p + runTimeG4_p + runTimeG5_p + runTimeG6_p + runTimeG7_p;
693 0 : T *tt=(T *)gridStore;
694 0 : grid.putStorage(tt,gDummy);
695 0 : }
696 : //
697 : //-----------------------------------------------------------------------------------
698 : // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding)
699 : //
700 0 : void AWVisResampler::GridToData(VBStore& vbs, const Array<Complex>& grid)
701 : {
702 : Int nDataChan, nDataPol, nGridPol, nGridChan, nx, ny,nw;//, nCFFreq;
703 : Int achan, apol, rbeg, rend;//, PolnPlane, ConjPlane;
704 0 : Vector<Float> sampling(2);//scaledSampling(2);
705 0 : Vector<Int> support(2),loc(3), iloc(4),tiloc(4);// scaledSupport(2);
706 0 : Vector<Double> pos(3), off(3);
707 :
708 0 : IPosition grdpos(4);
709 :
710 0 : Vector<Complex> norm(4,0.0);
711 0 : Complex phasor, nvalue;
712 0 : Vector<Int> cfShape=vbRow2CFBMap_p(0)->getStorage()(0,0,0)->getStorage()->shape().asVector();
713 0 : Vector<Double> pointingOffset((*vbRow2CFBMap_p(0)).getPointingOffset());
714 :
715 : // Vector<Int> convOrigin = (cfShape-1)/2;
716 0 : Vector<Int> convOrigin = (cfShape)/2;
717 0 : Double sinDPA=0.0, cosDPA=1.0, cfRefFreq;//cfScale=1.0
718 : // Int wndx = 0, fndx=0;
719 :
720 0 : rbeg=0;
721 0 : rend=vbs.nRow_p;
722 0 : rbeg = vbs.beginRow_p;
723 0 : rend = vbs.endRow_p;
724 0 : nx = grid.shape()[0]; ny = grid.shape()[1];
725 : //nw = cfShape[2];
726 0 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
727 :
728 0 : nDataPol = vbs.flagCube_p.shape()[0];
729 0 : nDataChan = vbs.flagCube_p.shape()[1];
730 :
731 : //
732 : // The following code reduces most array accesses to the simplest
733 : // possible to improve performance. However this made no
734 : // difference in the run-time performance compared to Vector,
735 : // Matrix and Cube indexing.
736 : //
737 : Bool Dummy;
738 0 : const Complex* __restrict__ gridStore = grid.getStorage(Dummy);
739 : (void)gridStore;
740 0 : Vector<Int> igrdpos(4);
741 0 : Double *freq=vbs.freq_p.getStorage(Dummy);
742 0 : Bool *rowFlag=vbs.rowFlag_p.getStorage(Dummy);
743 :
744 0 : Matrix<Double>& uvw=vbs.uvw_p;
745 0 : Cube<Complex>& visCube=vbs.visCube_p;
746 0 : Cube<Bool>& flagCube=vbs.flagCube_p;
747 :
748 0 : Vector<Int> gridInc, cfInc;
749 :
750 0 : cacheAxisIncrements(grid.shape().asVector(), gridInc_p);
751 : // cacheAxisIncrements(cfShape, cfInc_p);
752 : // Initialize the co-ordinates used for reading the CF values The
753 : // CFs are 4D arrays, with the last two axis degenerate (of length
754 : // 1). The last two axis were formerly the W-, and
755 : // Polarization-axis.
756 0 : iloc = 0;
757 : Bool finitePointingOffset=(
758 0 : (fabs(pointingOffset(0))>0) ||
759 0 : (fabs(pointingOffset(1))>0)
760 0 : );
761 0 : Int vbSpw = (vbs.vb_p)->spectralWindow();
762 :
763 0 : for(Int irow=rbeg; irow<rend; irow++) {
764 0 : if(!rowFlag[irow]) {
765 0 : CFBuffer& cfb = *vbRow2CFBMap_p(irow);
766 0 : Vector<Double> wVals, fVals; PolMapType mVals, mNdx, conjMVals, conjMNdx;
767 : Double fIncr, wIncr;
768 0 : cfb.getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
769 0 : nw = wVals.nelements();
770 : // nCFFreq = fVals.nelements()-1;
771 :
772 0 : for (Int ichan=0; ichan < nDataChan; ichan++) {
773 0 : achan=chanMap_p[ichan];
774 :
775 0 : if((achan>=0) && (achan<nGridChan)) {
776 : // lambda = C::c/freq[ichan];
777 0 : Double dataWVal = (vbs.vb_p->uvw()(irow)(2));
778 0 : Int wndx = cfb.nearestWNdx(abs(dataWVal)*freq[ichan]/C::c);
779 : //Int fndx = cfb.nearestFreqNdx(freq[ichan]);
780 0 : Int fndx = cfb.nearestFreqNdx(vbSpw,ichan);
781 :
782 : // cerr << "DG: " << fndx << " " << wndx << " " << ichan << " " << vbSpw << " " << freq[ichan] << endl;
783 :
784 : // cerr << "Grid: " << ichan << " " << freq[ichan] << " " << fndx << endl;
785 :
786 : // if (nw > 1) wndx=SynthesisUtils::nint((dataWVal*freq[ichan]/C::c)/wIncr-1);
787 : // if (nCFFreq > 0) fndx = SynthesisUtils::nint((freq[ichan])/fIncr-1);
788 : Float s;
789 : // CoordinateSystem cs;
790 : // cfb.getParams(cs,s,support(0),support(1),0,wndx,0);
791 0 : cfb.getParams(cfRefFreq,s,support(0),support(1),fndx,wndx,0);
792 0 : sampling(0) = sampling(1) = SynthesisUtils::nint(s);
793 :
794 : //cfScale = cfRefFreq/freq[ichan];
795 :
796 : // sampling[0] = SynthesisUtils::nint(sampling[0]*cfScale);
797 : // sampling[1] = SynthesisUtils::nint(sampling[1]*cfScale);
798 : // support[0] = SynthesisUtils::nint(support[0]/cfScale);
799 : // support[1] = SynthesisUtils::nint(support[1]/cfScale);
800 :
801 0 : sgrid(pos,loc,off,phasor,irow,uvw,dphase_p[irow],freq[ichan],
802 0 : uvwScale_p,offset_p,sampling);
803 :
804 : // iloc[2]=max(0, min(nw, loc[2]));
805 :
806 : Bool isOnGrid;
807 : //if ((isOnGrid=onGrid(nx, ny, nw, loc, support)))
808 : {
809 0 : for(Int ipol=0; ipol < nDataPol; ipol++) {
810 :
811 0 : if(!flagCube(ipol,ichan,irow)) {
812 0 : apol=polMap_p[ipol];
813 :
814 0 : if((apol>=0) && (apol<nGridPol)) {
815 0 : igrdpos[2]=apol; igrdpos[3]=achan;
816 0 : nvalue=0.0; norm(ipol)=0.0;
817 :
818 : // With VBRow2CFMap in use, CF for each pol. plane is a separate 2D Array.
819 0 : for (uInt mCol=0; mCol<conjMNdx[ipol].nelements(); mCol++)
820 : {
821 : int visGridElement, muellerElement;
822 : //
823 : // Get the pointer to the storage for the CF indexed by the Freq, W-term and Mueller
824 : // Element.
825 : //
826 0 : Complex* convFuncV=NULL;
827 : try
828 : {
829 0 : convFuncV = getConvFunc_p(vbs.paQuant_p.getValue("rad"),
830 : cfShape, support, muellerElement,
831 : cfb, dataWVal, fndx, wndx,
832 : //mNdx,conjMNdx,
833 : conjMNdx, mNdx,
834 : ipol, mCol);
835 : }
836 0 : catch (SynthesisFTMachineError& x)
837 : {
838 0 : LogIO log_l(LogOrigin("AWVisResampler[R&D]","GridToData"));
839 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
840 0 : }
841 : // Set the polarization plane of the gridded data to use for predicting with the CF from mCols column
842 0 : visGridElement=(int)(muellerElement%nDataPol);
843 0 : igrdpos[2]=polMap_p[visGridElement];
844 : //cerr << "DG: " << mCol << "-->" << visGridElement << "-->" << ipol << " " << polMap_p[ipol] << " " << polMap_p[visGridElement] << endl;
845 : //
846 : // Compute the incrmenets and center pixel for the current CF
847 : //
848 :
849 0 : if ((isOnGrid=onGrid(nx, ny, nw, loc, support))==false) break;
850 :
851 0 : cacheAxisIncrements(cfShape, cfInc_p);
852 0 : convOrigin = (cfShape)/2;
853 0 : if (finitePointingOffset)
854 0 : cachePhaseGrad_p(pointingOffset, cfShape, convOrigin, cfRefFreq, vbs.imRefFreq(),
855 0 : ((const Int)(vbs.vb_p)->spectralWindow()),((const Int)((vbs.vb_p)->fieldId())));
856 :
857 : // accumulateFromGrid() is a local C++ method with the inner loops. The include
858 : // file (FortanizedLoopsFromGrid.cc) has the interface code to call the inner
859 : // loops re-written in FORTRAN (in synthesis/fortran/faccumulateOnGrid.f)
860 :
861 : // accumulateFromGrid(nvalue, gridStore, igrdpos, convFuncV, dataWVal,
862 : // scaledSupport, scaledSampling, off, convOrigin,
863 : // cfShape, loc, phasor, sinDPA, cosDPA,
864 : // finitePointingOffset, cached_phaseGrad_p);
865 : #include <synthesis/TransformMachines/FortranizedLoopsFromGrid.inc>
866 :
867 : }
868 : // Zero divided by zero is NaN (IEEE standard)
869 0 : if (norm[ipol] != Complex(0.0)) visCube(ipol,ichan,irow)=nvalue/norm[ipol]; // Goes with FortranizedLoopsFromGrid.cc
870 : // visCube(ipol,ichan,irow)=nvalue; // Goes with FortranizedLoopsFromGrid.cc
871 : //if (casa::isNaN(nvalue))
872 : // {
873 : // cout << ipol << "," << ichan << "," << irow << "," << nvalue << "," << nDataChan << "," << nGridChan << "," << achan << endl;
874 : // //exit(0);
875 : // }
876 :
877 : //visCube(ipol,ichan,irow)=nvalue*conj(phasor)/norm(apol); // Goes with C++ loops
878 : // cerr << ipol << " " << ichan << " " << irow << " " << nvalue << " " << norm(apol) << " " << pointingOffset
879 : // << " " << qualifier_p << " " << ttt << " " << scaledSupport << endl;
880 : }
881 : }
882 : }
883 : }
884 : }
885 : }
886 : // junk++;
887 0 : }
888 : } // End row-loop
889 : // cerr << endl;
890 : // if (junk==20) exit(0);
891 0 : }
892 : //
893 : //-----------------------------------------------------------------------------------
894 : //
895 0 : void AWVisResampler::sgrid(Vector<Double>& pos, Vector<Int>& loc,
896 : Vector<Double>& off, Complex& phasor,
897 : const Int& irow, const Matrix<Double>& uvw,
898 : const Double& dphase, const Double& freq,
899 : const Vector<Double>& scale,
900 : const Vector<Double>& offset,
901 : const Vector<Float>& sampling)
902 : {
903 : Double phase;
904 0 : Vector<Double> uvw_l(3,0); // This allows gridding of weights
905 : // centered on the uv-origin
906 0 : if (uvw.nelements() > 0) for(Int i=0;i<3;i++) uvw_l[i]=uvw(i,irow);
907 :
908 0 : pos(2)=sqrt(abs(scale[2]*uvw_l(2)*freq/C::c))+offset[2];
909 0 : loc(2)=SynthesisUtils::nint(pos[2]);
910 0 : off(2)=0;
911 :
912 0 : for(Int idim=0;idim<2;idim++)
913 : {
914 0 : pos[idim]=scale[idim]*uvw_l(idim)*freq/C::c+(offset[idim]);
915 0 : loc[idim]=SynthesisUtils::nint(pos[idim]);
916 : // off[idim]=SynthesisUtils::nint((loc[idim]-pos[idim])*sampling[idim]+1);
917 0 : off[idim]=SynthesisUtils::nint((loc[idim]-pos[idim])*sampling[idim]);
918 : }
919 :
920 0 : if (dphase != 0.0)
921 : {
922 0 : phase=-2.0*C::pi*dphase*freq/C::c;
923 : Double sp,cp;
924 0 : SINCOS(phase,sp,cp);
925 : // phasor=Complex(cos(phase), sin(phase));
926 0 : phasor=Complex(cp,sp);
927 : }
928 : else
929 0 : phasor=Complex(1.0);
930 : // cerr << "### " << pos[0] << " " << offset[0] << " " << loc[0] << " " << off[0] << " " << uvw_l(0) << endl;
931 : // exit(0);
932 0 : }
933 : //
934 : //-----------------------------------------------------------------------------------
935 : //
936 0 : Bool AWVisResampler::reindex(const Vector<Int>& in, Vector<Int>& out,
937 : const Double& sinDPA, const Double& cosDPA,
938 : const Vector<Int>& Origin, const Vector<Int>& size)
939 : {
940 :
941 0 : Bool onGrid=false;
942 0 : Int ix=in[0], iy=in[1];
943 0 : if (sinDPA != 0.0)
944 : {
945 0 : ix = SynthesisUtils::nint(cosDPA*in[0] + sinDPA*in[1]);
946 0 : iy = SynthesisUtils::nint(-sinDPA*in[0] + cosDPA*in[1]);
947 : }
948 0 : out[0]=ix+Origin[0];
949 0 : out[1]=iy+Origin[1];
950 :
951 0 : onGrid = ((out[0] >= 0) && (out[0] < size[0]) &&
952 0 : (out[1] >= 0) && (out[1] < size[1]));
953 0 : if (!onGrid)
954 0 : cerr << "CF index out of range: " << out << " " << size << endl;
955 0 : return onGrid;
956 : }
957 :
958 :
959 : // void lineCFArea(const Int& th,
960 : // const Double& sinDPA,
961 : // const Double& cosDPA,
962 : // const Complex*__restrict__& convFuncV,
963 : // const Vector<Int>& cfShape,
964 : // const Vector<Int>& convOrigin,
965 : // const Int& cfInc,
966 : // Vector<Int>& iloc,
967 : // Vector<Int>& tiloc,
968 : // const Int* supportPtr,
969 : // const Float* samplingPtr,
970 : // const Double* offPtr,
971 : // Complex *cfAreaArrPtr)
972 : // {
973 : // cfAreaArrPtr[th]=0.0;
974 : // for(Int ix=-supportPtr[0]; ix <= supportPtr[0]; ix++)
975 : // {
976 : // iloc[0]=(Int)((samplingPtr[0]*ix+offPtr[0])-1);//+convOrigin[0];
977 : // tiloc=iloc;
978 : // if (reindex(iloc,tiloc,sinDPA, cosDPA,
979 : // convOrigin, cfShape))
980 : // {
981 : // wt = getFrom4DArray((const Complex * __restrict__ &)convFuncV,
982 : // tiloc,cfInc);
983 : // if (dataWVal > 0.0) wt = conj(wt);
984 : // cfAreaArrPtr[th] += wt;
985 : // }
986 : // }
987 : // }
988 :
989 0 : Complex AWVisResampler::getCFArea(Complex* __restrict__& convFuncV,
990 : Double& wVal,
991 : Vector<Int>& scaledSupport,
992 : Vector<Float>& scaledSampling,
993 : Vector<Double>& off,
994 : Vector<Int>& convOrigin,
995 : Vector<Int>& cfShape,
996 : Double& sinDPA,
997 : Double& cosDPA)
998 : {
999 0 : Vector<Int> iloc(4,0),tiloc(4);
1000 0 : Complex cfArea=0, wt;
1001 : Bool dummy;
1002 0 : Int *supportPtr=scaledSupport.getStorage(dummy);
1003 0 : Double *offPtr=off.getStorage(dummy);
1004 0 : Float *samplingPtr=scaledSampling.getStorage(dummy);
1005 0 : Int Nth=1;
1006 0 : Vector<Complex> cfAreaArr(Nth);
1007 0 : Complex *cfAreaArrPtr=cfAreaArr.getStorage(dummy);
1008 :
1009 0 : for(Int iy=-supportPtr[1]; iy <= supportPtr[1]; iy++)
1010 : {
1011 0 : iloc(1)=(Int)((samplingPtr[1]*iy+offPtr[1])-1);//+convOrigin[1];
1012 0 : for (Int th=0;th<Nth;th++)
1013 : {
1014 0 : cfAreaArr[th]=0.0;
1015 0 : for(Int ix=-supportPtr[0]; ix <= supportPtr[0]; ix++)
1016 : {
1017 0 : iloc[0]=(Int)((samplingPtr[0]*ix+offPtr[0])-1);//+convOrigin[0];
1018 0 : tiloc=iloc;
1019 0 : if (reindex(iloc,tiloc,sinDPA, cosDPA,
1020 : convOrigin, cfShape))
1021 : {
1022 0 : wt = getFrom4DArray((const Complex * __restrict__ &)convFuncV,
1023 0 : tiloc,cfInc_p);
1024 0 : if (wVal > 0.0) wt = conj(wt);
1025 0 : cfAreaArrPtr[th] += wt;
1026 : }
1027 : }
1028 : }
1029 0 : cfArea += sum(cfAreaArr);
1030 : }
1031 : // cerr << "cfArea: " << scaledSupport << " " << scaledSampling << " " << cfShape << " " << convOrigin << " " << cfArea << endl;
1032 0 : return cfArea;
1033 0 : }
1034 : using namespace casacore;
1035 : };// end namespace casa
|
