Line data Source code
1 : // -*- C++ -*-
2 : //# VisibilityResampler.cc: Implementation of the VisibilityResampler 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/VisibilityResampler.h>
31 : #include <synthesis/TransformMachines/Utils.h>
32 : #include <stdcasa/thread/AsynchronousTools.h>
33 : #include <fstream>
34 :
35 : using namespace casacore;
36 : namespace casa{
37 : #define NEED_UNDERSCORES
38 : #if defined(NEED_UNDERSCORES)
39 : #define gridengine gridengine_
40 : #define addtogrid addtogrid_
41 : #endif
42 : extern "C"{
43 : void gridengine(int *gnx,int *gny,int *gnpol,int *gnchan,
44 : double *grid,
45 : int *ncf, double *convfunc,
46 : float *sampling, int* support,
47 : int *loc, int* off, int *achan, int *apol,
48 : double * norm,
49 : float *phasor,
50 : int* ipol,int* ichan,int* irow,
51 : float* imgWts,int *dopsf,
52 : float* visCube,
53 : int *visCubePol, int *visCubeChan, int *visCubeRow);
54 : void addtogrid(double* grid, int* pos, float *val, double* wt,
55 : int *nx, int *ny, int *npol, int* nchan);
56 : }
57 : //
58 : //-----------------------------------------------------------------------------------
59 : //
60 : // VisibilityResampler& VisibilityResampler::operator=(const VisibilityResampler& other)
61 : // {
62 : // SynthesisUtils::SETVEC(uvwScale_p, other.uvwScale_p);
63 : // SynthesisUtils::SETVEC(offset_p, other.offset_p);
64 : // SynthesisUtils::SETVEC(dphase_p, other.dphase_p);
65 : // SynthesisUtils::SETVEC(chanMap_p, other.chanMap_p);
66 : // SynthesisUtils::SETVEC(polMap_p, other.polMap_p);
67 :
68 : // convFuncStore_p = other.convFuncStore_p;
69 : // myMutex_p = other.myMutex_p;
70 : // return *this;
71 : // }
72 : //
73 : //-----------------------------------------------------------------------------------
74 : // Re-sample the griddedData on the VisBuffer (a.k.a gridding)
75 : //
76 : // Template instantiations for re-sampling onto a double precision
77 : // or single precision grid.
78 : //
79 : template
80 : void VisibilityResampler::DataToGridImpl_p(Array<DComplex>& grid, VBStore& vbs,
81 : const Bool& dopsf, Matrix<Double>& sumwt,
82 : Bool /*useConjFreqCF*/); // __restrict__;
83 : template
84 : void VisibilityResampler::DataToGridImpl_p(Array<Complex>& grid, VBStore& vbs,
85 : const Bool& dopsf, Matrix<Double>& sumwt,
86 : Bool /*useConjFreqCF*/); // __restrict__;
87 :
88 : // template void VisibilityResampler::addTo4DArray(DComplex* store,const Int* iPos, Complex& val, Double& wt) __restrict__;
89 : // template void VisibilityResampler::addTo4DArray(Complex* store,const Int* iPos, Complex& val, Double& wt) __restrict__;
90 : //
91 : //-----------------------------------------------------------------------------------
92 : // Template implementation for DataToGrid
93 : //
94 : /*
95 : template <class T>
96 : void VisibilityResampler::DataToGridImpl_p(Array<T>& grid, VBStore& vbs, const Bool& dopsf,
97 : Matrix<Double>& sumwt)
98 : {
99 : Int nDataChan, nDataPol, nGridPol, nGridChan, nx, ny;
100 : Int achan, apol, rbeg, rend;
101 : Vector<Float> sampling(2);
102 : Vector<Int> support(2),loc(2), off(2), iloc(2);
103 : Vector<Double> pos(2);
104 :
105 : // IPosition grdpos(4);
106 : Vector<Int> igrdpos(4), gridIncrements;
107 :
108 : Double norm=0, wt, imgWt;
109 : Complex phasor, nvalue;
110 :
111 : rbeg = vbs.beginRow_p;
112 : rend = vbs.endRow_p;
113 : // cerr << rbeg << " " << rend << " " << vbs.nRow() << endl;
114 : nx = grid.shape()[0]; ny = grid.shape()[1];
115 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
116 : gridIncrements = grid.shape().asVector();
117 : nDataPol = vbs.flagCube_p.shape()[0];
118 : nDataChan = vbs.flagCube_p.shape()[1];
119 :
120 : sampling[0] = sampling[1] = convFuncStore_p.sampling[0];
121 : support(0) = convFuncStore_p.xSupport[0];
122 : support(1) = convFuncStore_p.ySupport[0];
123 :
124 : Bool Dummy, gDummy;
125 : T __restrict__ *gridStore = grid.getStorage(gDummy);
126 : Int * __restrict__ iPosPtr = igrdpos.getStorage(Dummy);
127 : const Int * __restrict__ iPosPtrConst = iPosPtr;
128 : Double *__restrict__ convFunc=(*(convFuncStore_p.rdata)).getStorage(Dummy);
129 : Double * __restrict__ freq=vbs.freq_p.getStorage(Dummy);
130 : Bool * __restrict__ rowFlag=vbs.rowFlag_p.getStorage(Dummy);
131 :
132 : Float * __restrict__ imagingWeight = vbs.imagingWeight_p.getStorage(Dummy);
133 : Double * __restrict__ uvw = vbs.uvw_p.getStorage(Dummy);
134 : Bool * __restrict__ flagCube = vbs.flagCube_p.getStorage(Dummy);
135 : Complex * __restrict__ visCube = vbs.visCube_p.getStorage(Dummy);
136 : Double * __restrict__ scale = uvwScale_p.getStorage(Dummy);
137 : Double * __restrict__ offset = offset_p.getStorage(Dummy);
138 : Float * __restrict__ samplingPtr = sampling.getStorage(Dummy);
139 : Double * __restrict__ posPtr=pos.getStorage(Dummy);
140 : Int * __restrict__ locPtr=loc.getStorage(Dummy);
141 : Int * __restrict__ offPtr=off.getStorage(Dummy);
142 : Double * __restrict__ sumwtPtr = sumwt.getStorage(Dummy);
143 : Int * __restrict__ ilocPtr=iloc.getStorage(Dummy);
144 : Int * __restrict__ supportPtr = support.getStorage(Dummy);
145 : Int nDim = vbs.uvw_p.shape()[0];
146 :
147 : // cacheAxisIncrements(nx,ny,nGridPol, nGridChan);
148 : cacheAxisIncrements(grid.shape().asVector());
149 :
150 : for(Int irow=rbeg; irow < rend; irow++){ // For all rows
151 :
152 : if(!rowFlag[irow]){ // If the row is not flagged
153 :
154 : for(Int ichan=0; ichan< nDataChan; ichan++){ // For all channels
155 :
156 : // if (vbs.imagingWeight(ichan,irow)!=0.0) { // If weights are not zero
157 : if (imagingWeight[ichan+irow*nDataChan]!=0.0) { // If weights are not zero
158 : achan=chanMap_p(ichan);
159 :
160 : if((achan>=0) && (achan<nGridChan)) { // If selected channels are valid
161 :
162 : // sgrid(pos,loc,off, phasor, irow,
163 : // vbs.uvw,dphase_p[irow], vbs.freq[ichan],
164 : // uvwScale_p, offset_p, sampling);
165 : sgrid(nDim,posPtr,locPtr,offPtr, phasor, irow,
166 : uvw,dphase_p[irow], freq[ichan],
167 : scale, offset, samplingPtr);
168 :
169 : if (onGrid(nx, ny, loc, support)) { // If the data co-ords. are with-in the grid
170 :
171 : for(Int ipol=0; ipol< nDataPol; ipol++) { // For all polarizations
172 : // if((!vbs.flagCube(ipol,ichan,irow))){ // If the pol. & chan. specific
173 : if((!flagCube[ipol+ichan*nDataPol+irow*nDataChan*nDataPol])){
174 : apol=polMap_p(ipol);
175 : if ((apol>=0) && (apol<nGridPol)) {
176 : // igrdpos(2)=apol; igrdpos(3)=achan;
177 : iPosPtr[2]=apol; iPosPtr[3]=achan;
178 : norm=0.0;
179 :
180 : imgWt=imagingWeight[ichan+irow*nDataChan];
181 :
182 : Int idopsf = (dopsf==true);
183 : Int ncf=(*convFuncStore_p.rdata).shape()(0);
184 : Int nrow = vbs.nRow();
185 : gridengine(&nx, &ny, &nGridPol, &nGridChan, (double*)gridStore,
186 : &ncf, (double*)convFunc, samplingPtr, supportPtr,
187 : locPtr, offPtr, &achan, &apol, &norm, (float *)&phasor,
188 : &ipol, &ichan, &irow, imagingWeight, &idopsf,
189 : (float *)visCube, &nDataPol, &nDataChan,&nrow);
190 :
191 : // if(dopsf) nvalue=Complex(imgWt);
192 : // else nvalue=imgWt*
193 : // // (vbs.visCube(ipol,ichan,irow)*phasor);
194 : // (visCube[ipol+ichan*nDataPol+irow*nDataPol*nDataChan]*phasor);
195 :
196 : // for(Int iy=-supportPtr[1]; iy <= supportPtr[1]; iy++)
197 : // {
198 : // ilocPtr[1]=abs((int)(samplingPtr[1]*iy+offPtr[1]));
199 : // // igrdpos(1)=loc(1)+iy;
200 : // iPosPtr[1]=locPtr[1]+iy;
201 : // // wt = convFunc[ilocPtr[1]];
202 : // for(Int ix=-supportPtr[0]; ix <= supportPtr[0]; ix++)
203 : // {
204 : // ilocPtr[0]=abs((int)(samplingPtr[0]*ix+offPtr[0]));
205 : // wt = convFunc[iloc[0]]*convFunc[iloc[1]];
206 : // // wt *= convFunc[ilocPtr[0]];
207 :
208 : // //igrdpos(0)=loc(0)+ix;
209 : // iPosPtr[0]=locPtr[0]+ix;
210 :
211 :
212 : // // grid(grdpos) += nvalue*wt;
213 :
214 : // // gridStore[iPosPtr[0] +
215 : // // iPosPtr[1]*incPtr_p[1] +
216 : // // iPosPtr[2]*incPtr_p[2] +
217 : // // iPosPtr[3]*incPtr_p[3]].real() += (nvalue.real()*wt);
218 : // // gridStore[iPosPtr[0] +
219 : // // iPosPtr[1]*incPtr_p[1] +
220 : // // iPosPtr[2]*incPtr_p[2] +
221 : // // iPosPtr[3]*incPtr_p[3]].imag() += (nvalue.imag()*wt);
222 :
223 : // // The following uses raw index on the 4D grid
224 : // // addTo4DArray(gridStore,iPosPtr,nvalue,wt);
225 :
226 : // addtogrid((double *)gridStore, iPosPtr, (float *)&nvalue, (double *)&wt,
227 : // &nx, &ny, &nGridPol, &nGridChan);
228 :
229 : // norm+=wt;
230 : // }
231 : // }
232 : sumwtPtr[apol+achan*nGridPol]+=imgWt*norm;
233 : // sumwt(apol,achan)+=imgWt*norm;
234 : }
235 : }
236 : }
237 : }
238 : }
239 : }
240 : }
241 : }
242 : }
243 : // T *tt=(T *)gridStore;
244 : // grid.putStorage(tt,gDummy);
245 : }
246 : */
247 : //
248 : //-----------------------------------------------------------------------------------
249 : // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding)
250 : //
251 0 : void VisibilityResampler::GridToData(VBStore& vbs, const Array<Complex>& grid)
252 : // void VisibilityResampler::GridToData(VBStore& vbs, Array<Complex>& grid)
253 : {
254 : using casacore::operator*;
255 :
256 : Int nDataChan, nDataPol, nGridPol, nGridChan, nx, ny;
257 : Int achan, apol, rbeg, rend;
258 0 : Vector<Float> sampling(2);
259 0 : Vector<Int> support(2),loc(2), off(2), iloc(2);
260 0 : Vector<Double> pos(2);
261 :
262 0 : IPosition grdpos(4);
263 :
264 0 : Double norm=0, wt;
265 0 : Complex phasor, nvalue;
266 :
267 : // rbeg=0;
268 : // rend=vbs.nRow_p;
269 0 : rbeg = vbs.beginRow_p;
270 0 : rend = vbs.endRow_p;
271 : // cerr << rbeg << " " << rend << " " << vbs.nRow() << endl;
272 0 : nx = grid.shape()[0]; ny = grid.shape()[1];
273 0 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
274 :
275 0 : nDataPol = vbs.flagCube_p.shape()[0];
276 0 : nDataChan = vbs.flagCube_p.shape()[1];
277 :
278 0 : sampling[0] = sampling[1] = convFuncStore_p.sampling[0];
279 0 : support(0) = convFuncStore_p.xSupport[0];
280 0 : support(1) = convFuncStore_p.ySupport[0];
281 :
282 : Bool Dummy,vbcDummy;
283 0 : const Complex *__restrict__ gridStore = grid.getStorage(Dummy);
284 0 : Vector<Int> igrdpos(4);
285 0 : Int *__restrict__ iPosPtr = igrdpos.getStorage(Dummy);
286 0 : const Int *__restrict__ iPosPtrConst = iPosPtr;
287 0 : Double *__restrict__ convFunc=(*(convFuncStore_p.rdata)).getStorage(Dummy);
288 0 : Double *__restrict__ freq=vbs.freq_p.getStorage(Dummy);
289 0 : Bool *__restrict__ rowFlag=vbs.rowFlag_p.getStorage(Dummy);
290 :
291 : //UNUSED: Float * __restrict__ imagingWeight = vbs.imagingWeight_p.getStorage(Dummy);
292 0 : Double * __restrict__ uvw = vbs.uvw_p.getStorage(Dummy);
293 0 : Bool * __restrict__ flagCube = vbs.flagCube_p.getStorage(Dummy);
294 0 : Complex * __restrict__ visCube = vbs.visCube_p.getStorage(vbcDummy);
295 0 : Double * __restrict__ scale = uvwScale_p.getStorage(Dummy);
296 0 : Double * __restrict__ offset = offset_p.getStorage(Dummy);
297 0 : Int * __restrict__ supportPtr = support.getStorage(Dummy);
298 0 : Float * __restrict__ samplingPtr = sampling.getStorage(Dummy);
299 0 : Double * __restrict__ posPtr=pos.getStorage(Dummy);
300 0 : Int * __restrict__ locPtr=loc.getStorage(Dummy);
301 0 : Int * __restrict__ offPtr=off.getStorage(Dummy);
302 0 : Int * __restrict__ ilocPtr=iloc.getStorage(Dummy);
303 0 : Int nDim = vbs.uvw_p.shape()(0);
304 :
305 : // cacheAxisIncrements(nx,ny,nGridPol, nGridChan);
306 0 : cacheAxisIncrements(grid.shape().asVector());
307 : //UNUSED: Int xSupport_l, ySupport_l;
308 : //UNUSED: Float sampling_l;
309 0 : for(Int irow=rbeg; irow < rend; irow++) {
310 0 : if(!rowFlag[irow]) {
311 0 : for (Int ichan=0; ichan < nDataChan; ichan++) {
312 0 : achan=chanMap_p(ichan);
313 :
314 0 : if((achan>=0) && (achan<nGridChan)) {
315 : // sgrid(pos,loc,off,phasor,irow,vbs.uvw,
316 : // dphase_p[irow],vbs.freq[ichan],
317 : // uvwScale_p,offset_p,sampling);
318 :
319 0 : sgrid(nDim,posPtr,locPtr,offPtr,phasor,irow,uvw,
320 0 : dphase_p[irow],freq[ichan],
321 : scale,offset,samplingPtr);
322 :
323 0 : if (onGrid(nx, ny, loc, support)) {
324 0 : for(Int ipol=0; ipol < nDataPol; ipol++) {
325 :
326 0 : if(!flagCube[ipol+ichan*nDataPol+irow*nDataChan*nDataPol]) {
327 0 : apol=polMap_p(ipol);
328 :
329 0 : if((apol>=0) && (apol<nGridPol)) {
330 : // igrdpos(2)=apol; igrdpos(3)=achan;
331 0 : iPosPtr[2]=apol; iPosPtr[3]=achan;
332 0 : nvalue=0.0;
333 0 : norm=0.0;
334 :
335 0 : for(Int iy=-supportPtr[1]; iy <= supportPtr[1]; iy++)
336 : {
337 0 : ilocPtr[1]=abs((Int)(samplingPtr[1]*iy+offPtr[1]));
338 : // igrdpos(1)=loc(1)+iy;
339 0 : iPosPtr[1]=locPtr[1]+iy;
340 : // wt = convFunc[ilocPtr[1]];
341 0 : for(Int ix=-supportPtr[0]; ix <= supportPtr[0]; ix++)
342 : {
343 0 : ilocPtr[0]=abs((Int)(samplingPtr[0]*ix+offPtr[0]));
344 : // igrdpos(0)=loc(0)+ix;
345 0 : iPosPtr[0]=locPtr[0]+ix;
346 0 : wt=convFunc[ilocPtr[0]]*convFunc[ilocPtr[1]];
347 : // wt *= convFunc[ilocPtr[0]];
348 0 : norm+=wt;
349 : // nvalue+=wt*grid(grdpos);
350 : // The following uses raw index on the 4D grid
351 0 : nvalue+=wt*getFrom4DArray(gridStore,iPosPtrConst);
352 : }
353 : }
354 0 : visCube[ipol+ichan*nDataPol+irow*nDataChan*nDataPol]=(nvalue*conj(phasor))/norm;
355 : }
356 : }
357 : }
358 : }
359 : }
360 : }
361 : }
362 : }
363 0 : Complex *tt=(Complex *) visCube;
364 0 : vbs.visCube_p.putStorage(tt,vbcDummy);
365 0 : }
366 : //
367 : //-----------------------------------------------------------------------------------
368 : //
369 : // void VisibilityResampler::sgrid(Int& uvwDim,Double* __restrict__ pos,
370 : // Int* __restrict__ loc,
371 : // Int* __restrict__ off,
372 : // Complex& phasor, const Int& irow,
373 : // // const Matrix<Double>& __restrict__ uvw,
374 : // const Double* __restrict__ uvw,
375 : // const Double& __restrict__ dphase,
376 : // const Double& __restrict__ freq,
377 : // const Double* __restrict__ scale,
378 : // const Double* __restrict__ offset,
379 : // const Float* __restrict__ sampling) __restrict__
380 : // // const Vector<Double>& __restrict__ scale,
381 : // // const Vector<Double>& __restrict__ offset,
382 : // // const Vector<Float>& __restrict__ sampling) __restrict__
383 : // {
384 : // Double phase;
385 : // // Int ndim=pos.shape()(0);
386 : // Int ndim=2;
387 :
388 : // for(Int idim=0;idim<ndim;idim++)
389 : // {
390 : // pos[idim]=scale[idim]*uvw[idim+irow*uvwDim]*freq/C::c+offset[idim];
391 : // loc[idim]=(Int)std::floor(pos[idim]+0.5);
392 : // off[idim]=(Int)std::floor(((loc[idim]-pos[idim])*sampling[idim])+0.5);
393 : // }
394 :
395 : // if (dphase != 0.0)
396 : // {
397 : // phase=-2.0*C::pi*dphase*freq/C::c;
398 : // phasor=Complex(cos(phase), sin(phase));
399 : // }
400 : // else
401 : // phasor=Complex(1.0);
402 : // }
403 0 : void VisibilityResampler::ComputeResiduals(VBStore& vbs)
404 : {
405 0 : Int rbeg = vbs.beginRow_p, rend = vbs.endRow_p;
406 0 : IPosition vbDataShape=vbs.modelCube_p.shape();
407 0 : IPosition start(vbDataShape), last(vbDataShape);
408 0 : start=0; start(2)=rbeg;
409 0 : last(2)=rend; //last=last-1;
410 :
411 0 : if (!vbs.useCorrected_p)
412 : {
413 0 : for(uInt ichan = start(0); ichan < last(0); ichan++)
414 0 : for(uInt ipol = start(1); ipol < last(1); ipol++)
415 0 : for(uInt irow = start(2); irow < last(2); irow++)
416 0 : vbs.modelCube_p(ichan,ipol,irow) = vbs.modelCube_p(ichan,ipol,irow) - vbs.visCube_p(ichan,ipol,irow);
417 : // vbs.modelCube_p(ichan,ipol,irow) = vbs.visCube_p(ichan,ipol,irow) - vbs.modelCube_p(ichan,ipol,irow);
418 : }
419 : else
420 : {
421 0 : for(uInt ichan = start(0); ichan < last(0); ichan++)
422 0 : for(uInt ipol = start(1); ipol < last(1); ipol++)
423 0 : for(uInt irow = start(2); irow < last(2); irow++)
424 0 : vbs.modelCube_p(ichan,ipol,irow) = vbs.modelCube_p(ichan,ipol,irow) - vbs.correctedCube_p(ichan,ipol,irow);
425 : //vbs.modelCube_p(ichan,ipol,irow) = vbs.correctedCube_p(ichan,ipol,irow) - vbs.modelCube_p(ichan,ipol,irow);
426 : }
427 :
428 : // Slicer mySlice(start,last,Slicer::endIsLast);
429 : // Cube<Complex> slicedModel, slicedData, slicedCorrected;
430 : // if (!vbs.useCorrected_p)
431 : // {
432 : // {
433 : // async::MutexLocker tt(*myMutex_p);
434 : // slicedModel = Cube<Complex>(vbs.modelCube_p(mySlice));
435 : // slicedData = Cube<Complex>(vbs.visCube_p(mySlice));
436 : // }
437 : // slicedModel -= slicedData;
438 : // }
439 : // else
440 : // {
441 : // {
442 : // async::MutexLocker tt(*myMutex_p);
443 : // slicedModel = Cube<Complex>(vbs.modelCube_p(mySlice));
444 : // slicedCorrected = Cube<Complex>(vbs.correctedCube_p(mySlice));
445 : // }
446 : // slicedModel -= slicedCorrected;
447 : // }
448 0 : }
449 :
450 :
451 :
452 :
453 : #define DoOld 1
454 :
455 : template <class T>
456 0 : void VisibilityResampler::DataToGridImpl_p(Array<T>& grid, VBStore& vbs, const Bool& dopsf,
457 : Matrix<Double>& sumwt,
458 : Bool /*useConjFreqCF*/)
459 : {
460 : using casacore::operator*;
461 :
462 : static Bool beenThereDoneThat = false;
463 0 : if (! beenThereDoneThat){
464 : #if DoOld
465 0 : cerr << "==> Doing it the old way" << endl;
466 : #else
467 : cerr << "==> Doing it the new way" << endl;
468 : #endif
469 0 : beenThereDoneThat = true;
470 : }
471 :
472 : Int nDataChan, nDataPol, nGridPol, nGridChan, nx, ny;
473 : Int achan, apol, rbeg, rend;
474 0 : Vector<Float> sampling(2);
475 0 : Vector<Int> support(2),loc(2), off(2), iloc(2);
476 0 : Vector<Double> pos(2);
477 :
478 : // IPosition grdpos(4);
479 0 : Vector<Int> igrdpos(4);
480 :
481 0 : Double norm=0, wt, imgWt;
482 0 : Complex phasor, nvalue;
483 :
484 0 : rbeg = vbs.beginRow_p;
485 0 : rend = vbs.endRow_p;
486 : // cerr << rbeg << " " << rend << " " << vbs.nRow() << endl;
487 0 : nx = grid.shape()[0]; ny = grid.shape()[1];
488 0 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
489 :
490 0 : nDataPol = vbs.flagCube_p.shape()[0];
491 0 : nDataChan = vbs.flagCube_p.shape()[1];
492 :
493 0 : sampling[0] = sampling[1] = convFuncStore_p.sampling[0];
494 0 : support(0) = convFuncStore_p.xSupport[0];
495 0 : support(1) = convFuncStore_p.ySupport[0];
496 :
497 : Bool Dummy,gDummy;
498 0 : T * __restrict__ gridStore = grid.getStorage(gDummy);
499 0 : Int * __restrict__ iPosPtr = igrdpos.getStorage(Dummy);
500 0 : Double *__restrict__ convFunc=(*(convFuncStore_p.rdata)).getStorage(Dummy);
501 0 : Double * __restrict__ freq=vbs.freq_p.getStorage(Dummy);
502 0 : Bool * __restrict__ rowFlag=vbs.rowFlag_p.getStorage(Dummy);
503 :
504 0 : Float * __restrict__ imagingWeight = vbs.imagingWeight_p.getStorage(Dummy);
505 0 : Double * __restrict__ uvw = vbs.uvw_p.getStorage(Dummy);
506 0 : Bool * __restrict__ flagCube = vbs.flagCube_p.getStorage(Dummy);
507 0 : Complex * __restrict__ visCube = vbs.visCube_p.getStorage(Dummy);
508 0 : Double * __restrict__ scale = uvwScale_p.getStorage(Dummy);
509 0 : Double * __restrict__ offset = offset_p.getStorage(Dummy);
510 0 : Float * __restrict__ samplingPtr = sampling.getStorage(Dummy);
511 0 : Double * __restrict__ posPtr=pos.getStorage(Dummy);
512 0 : Int * __restrict__ locPtr=loc.getStorage(Dummy);
513 0 : Int * __restrict__ offPtr=off.getStorage(Dummy);
514 : //UNUSED: Double * __restrict__ sumwtPtr = sumwt.getStorage(Dummy);
515 0 : Int nDim = vbs.uvw_p.shape()[0];
516 :
517 : // cacheAxisIncrements(nx,ny,nGridPol, nGridChan);
518 0 : cacheAxisIncrements(grid.shape().asVector());
519 :
520 : // Vector<Int> gridIncrements = grid.shape().asVector();
521 0 : Vector<Int> gridShape = grid.shape().asVector();
522 0 : Vector<Int> gridIncrements (gridShape.nelements());
523 :
524 0 : gridIncrements[0] = 1;
525 0 : for (uint i = 1; i < gridShape.nelements(); i++){
526 0 : gridIncrements [i] = gridIncrements[i-1] * gridShape[i-1];
527 : }
528 0 : Vector<Double> convolutionLookupX (2 * support[0] + 1, 0.0);
529 0 : Vector<Double> convolutionLookupY (2 * support[1] + 1, 0.0);
530 : //UNUSED: const Double * const pConvolutionLookupY0 = convolutionLookupY.getStorage (Dummy);
531 : //UNUSED const Double * const pConvolutionLookupX0 = convolutionLookupX.getStorage (Dummy);
532 : //UNUSED: const Double * const pConvolutionLookupXEnd = pConvolutionLookupX0 + convolutionLookupX.size();
533 :
534 0 : for(Int irow=rbeg; irow < rend; irow++){ // For all rows
535 :
536 0 : if(!rowFlag[irow]){ // If the row is not flagged
537 :
538 0 : for(Int ichan=0; ichan< nDataChan; ichan++){ // For all channels
539 :
540 : // if (vbs.imagingWeight(ichan,irow)!=0.0) { // If weights are not zero
541 0 : if (imagingWeight[ichan+irow*nDataChan]!=0.0) { // If weights are not zero
542 0 : achan=chanMap_p[ichan];
543 :
544 0 : if((achan>=0) && (achan<nGridChan)) { // If selected channels are valid
545 :
546 : // sgrid(pos,loc,off, phasor, irow,
547 : // vbs.uvw,dphase_p[irow], vbs.freq[ichan],
548 : // uvwScale_p, offset_p, sampling);
549 0 : sgrid(nDim,posPtr,locPtr,offPtr, phasor, irow,
550 0 : uvw,dphase_p[irow], freq[ichan],
551 : scale, offset, samplingPtr);
552 :
553 0 : if (onGrid(nx, ny, loc, support)) { // If the data co-ords. are with-in the grid
554 :
555 0 : Double convolutionSumX = 0;
556 0 : for (int ix = - support [0], ii = 0; ix <= support [0]; ix ++, ii++){
557 0 : Int iConv = abs(int(sampling[0] * ix + off[0]));
558 0 : convolutionLookupX [ii] = convFunc[iConv];
559 0 : convolutionSumX += convFunc[iConv];
560 : }
561 :
562 0 : Double convolutionSumY= 0;
563 0 : for (int iy = - support [1], ii = 0; iy <= support [1]; iy ++, ii++){
564 0 : Int iConv = abs(int(sampling[1] * iy + off[1]));
565 0 : convolutionLookupY [ii] = convFunc[iConv];
566 0 : convolutionSumY += convFunc[iConv];
567 : }
568 :
569 : //UNUSED: Double Norm = convolutionSumX * convolutionSumY;
570 :
571 :
572 0 : for(Int ipol=0; ipol< nDataPol; ipol++) { // For all polarizations
573 : // if((!vbs.flagCube(ipol,ichan,irow))){ // If the pol. & chan. specific
574 0 : if((!flagCube[ipol+ichan*nDataPol+irow*nDataChan*nDataPol])){
575 0 : apol=polMap_p(ipol);
576 0 : if ((apol>=0) && (apol<nGridPol)) {
577 0 : igrdpos[2]=apol; igrdpos[3]=achan;
578 :
579 0 : norm=0.0;
580 :
581 0 : imgWt=imagingWeight[ichan+irow*nDataChan];
582 0 : if(dopsf) nvalue=Complex(imgWt);
583 0 : else nvalue=imgWt*
584 : // (vbs.visCube(ipol,ichan,irow)*phasor);
585 0 : (visCube[ipol+ichan*nDataPol+irow*nDataPol*nDataChan]*phasor);
586 :
587 : #if DoOld
588 :
589 : // Original Inner Loop
590 :
591 : //vector<const T*> oldAddresses;
592 :
593 : // off[idim]=(Int)std::floor(((loc[idim]-pos[idim])*sampling[idim])+0.5); un
594 :
595 0 : for(Int iy=-support[1]; iy <= support[1]; iy++)
596 : {
597 0 : iloc(1)=abs((int)(sampling[1]*iy+off[1]));
598 0 : igrdpos[1]=loc[1]+iy;
599 0 : for(Int ix=-support[0]; ix <= support[0]; ix++)
600 : {
601 0 : iloc[0]=abs((int)(sampling[0]*ix+off[0]));
602 0 : wt = convFunc[iloc[0]]*convFunc[iloc[1]];
603 :
604 0 : igrdpos[0]=loc[0]+ix;
605 : // grid(grdpos) += nvalue*wt;
606 :
607 : // The following uses raw index on the 4D grid
608 0 : addTo4DArray(gridStore,iPosPtr,nvalue,wt);
609 :
610 : //oldAddresses.push_back (& gridStore[igrdpos[0] + igrdpos[1]*gridIncrements[1] + igrdpos[2]*gridIncrements[2] +
611 : // igrdpos[3]*gridIncrements[3]]);
612 :
613 0 : norm+=wt;
614 : }
615 : }
616 : #else
617 :
618 : // New Inner Loop
619 :
620 : //vector<const T*> newAddresses;
621 :
622 :
623 : const Int X = 0;
624 : const Int Y = 1;
625 : const Int Z1 = 2; // channel
626 : const Int Z2 = 3; // polarization
627 :
628 : Int gridZ1 = igrdpos[Z1]; // Third grid coordinate, Z1
629 : Int gridZ2 = igrdpos[Z2]; // Fourth grid coordinate, Z2
630 :
631 : T * gridStoreZ1Z2 = gridStore + gridZ1 * gridIncrements [Z1] +
632 : gridZ2 * gridIncrements [Z2];
633 : // Position of origin of xy plane specified by Z1, Z2
634 :
635 : T * gridStoreYZ1Z2 =
636 : gridStoreZ1Z2 + gridIncrements [Y] * (loc[Y] - support[Y] - 1);
637 : // Position of origin of lower left corner of rectangle
638 : // of XY plane used in convolution
639 :
640 : //const Int offX = off[X];
641 : //const Int offY = off[Y];
642 : const Int yIncrement = gridIncrements [Y];
643 : const Int x0 = (loc[X] - support[X]) - 1;
644 : const Int xMax = support[X];
645 : const Int nX = xMax * 2 + 1;
646 : const Int yMax = support[Y];
647 : const Int nY = yMax * 2 + 1;
648 : const Double * pConvolutionLookupY = pConvolutionLookupY0;
649 :
650 : for(Int iy=0; iy < nY ; iy ++)
651 : {
652 : const Double convFuncY = * pConvolutionLookupY ++;
653 : const Double * __restrict__ pConvolutionLookupX = pConvolutionLookupX0;
654 :
655 : gridStoreYZ1Z2 += yIncrement;
656 : T * __restrict__ gridStoreXYZ1Z2 = gridStoreYZ1Z2 + x0;
657 :
658 : for (const Double * pConvolutionLookupX = pConvolutionLookupX0;
659 : pConvolutionLookupX != pConvolutionLookupXEnd;
660 : pConvolutionLookupX ++){
661 : // for(Int ix=0; ix < nX; ix ++)
662 : // {
663 : wt = (* pConvolutionLookupX) * convFuncY;
664 :
665 : gridStoreXYZ1Z2 += 1;
666 : //newAddresses.push_back (gridStoreXYZ1Z2);
667 : * gridStoreXYZ1Z2 += (nvalue * wt);
668 : //multiplyAndAdd (gridStoreXYZ1Z2, nvalue, wt);
669 :
670 : //norm += wt;
671 : }
672 : }
673 :
674 : norm = Norm;
675 :
676 :
677 : #endif
678 :
679 : // sumwtPtr[apol+achan*nGridPol]+=imgWt*norm;
680 0 : sumwt(apol,achan)+=imgWt*norm;
681 : }
682 : }
683 : }
684 : }
685 : }
686 : }
687 : }
688 : }
689 : }
690 0 : T *tt=(T *)gridStore;
691 0 : grid.putStorage(tt,gDummy);
692 0 : }
693 :
694 :
695 : using namespace casacore;
696 : };// end namespace casa
|
