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: 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 : #include <synthesis/TransformMachines/SynthesisError.h>
30 : #include <synthesis/TransformMachines2/VisibilityResampler.h>
31 : #include <synthesis/TransformMachines2/Utils.h>
32 : #include <stdcasa/thread/AsynchronousTools.h>
33 : #include <fstream>
34 :
35 : using namespace casacore;
36 : namespace casa{
37 : using namespace refim;
38 : #define NEED_UNDERSCORES
39 : #if defined(NEED_UNDERSCORES)
40 : #define gridengine gridengine_
41 : #define addtogrid addtogrid_
42 : #endif
43 : extern "C"{
44 : void gridengine(int *gnx,int *gny,int *gnpol,int *gnchan,
45 : double *grid,
46 : int *ncf, double *convfunc,
47 : float *sampling, int* support,
48 : int *loc, int* off, int *achan, int *apol,
49 : double * norm,
50 : float *phasor,
51 : int* ipol,int* ichan,int* irow,
52 : float* imgWts,int *dopsf,
53 : float* visCube,
54 : int *visCubePol, int *visCubeChan, int *visCubeRow);
55 : void addtogrid(double* grid, int* pos, float *val, double* wt,
56 : int *nx, int *ny, int *npol, int* nchan);
57 : }
58 : //
59 : //-----------------------------------------------------------------------------------
60 : //
61 : // VisibilityResampler& VisibilityResampler::operator=(const VisibilityResampler& other)
62 : // {
63 : // SynthesisUtils::SETVEC(uvwScale_p, other.uvwScale_p);
64 : // SynthesisUtils::SETVEC(offset_p, other.offset_p);
65 : // SynthesisUtils::SETVEC(dphase_p, other.dphase_p);
66 : // SynthesisUtils::SETVEC(chanMap_p, other.chanMap_p);
67 : // SynthesisUtils::SETVEC(polMap_p, other.polMap_p);
68 :
69 : // convFuncStore_p = other.convFuncStore_p;
70 : // myMutex_p = other.myMutex_p;
71 : // return *this;
72 : // }
73 : //
74 : //-----------------------------------------------------------------------------------
75 : // Re-sample the griddedData on the VisBuffer (a.k.a gridding)
76 : //
77 : // Template instantiations for re-sampling onto a double precision
78 : // or single precision grid.
79 : //
80 : template
81 : void VisibilityResampler::DataToGridImpl_p(Array<DComplex>& grid, VBStore& vbs,
82 : const Bool& dopsf, Matrix<Double>& sumwt,
83 : Bool /*useConjFreqCF*/); // __restrict__;
84 : template
85 : void VisibilityResampler::DataToGridImpl_p(Array<Complex>& grid, VBStore& vbs,
86 : const Bool& dopsf, Matrix<Double>& sumwt,
87 : Bool /*useConjFreqCF*/); // __restrict__;
88 :
89 : // template void VisibilityResampler::addTo4DArray(DComplex* store,const Int* iPos, Complex& val, Double& wt) __restrict__;
90 : // template void VisibilityResampler::addTo4DArray(Complex* store,const Int* iPos, Complex& val, Double& wt) __restrict__;
91 : //
92 : //-----------------------------------------------------------------------------------
93 : // Template implementation for DataToGrid
94 : //
95 : /*
96 : template <class T>
97 : void VisibilityResampler::DataToGridImpl_p(Array<T>& grid, VBStore& vbs, const Bool& dopsf,
98 : Matrix<Double>& sumwt)
99 : {
100 : Int nDataChan, nDataPol, nGridPol, nGridChan, nx, ny;
101 : Int achan, apol, rbeg, rend;
102 : Vector<Float> sampling(2);
103 : Vector<Int> support(2),loc(2), off(2), iloc(2);
104 : Vector<Double> pos(2);
105 :
106 : // IPosition grdpos(4);
107 : Vector<Int> igrdpos(4), gridIncrements;
108 :
109 : Double norm=0, wt, imgWt;
110 : Complex phasor, nvalue;
111 :
112 : rbeg = vbs.beginRow_p;
113 : rend = vbs.endRow_p;
114 : // cerr << rbeg << " " << rend << " " << vbs.nRow() << endl;
115 : nx = grid.shape()[0]; ny = grid.shape()[1];
116 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
117 : gridIncrements = grid.shape().asVector();
118 : nDataPol = vbs.flagCube_p.shape()[0];
119 : nDataChan = vbs.flagCube_p.shape()[1];
120 :
121 : sampling[0] = sampling[1] = convFuncStore_p.sampling[0];
122 : support(0) = convFuncStore_p.xSupport[0];
123 : support(1) = convFuncStore_p.ySupport[0];
124 :
125 : Bool Dummy, gDummy;
126 : T __restrict__ *gridStore = grid.getStorage(gDummy);
127 : Int * __restrict__ iPosPtr = igrdpos.getStorage(Dummy);
128 : const Int * __restrict__ iPosPtrConst = iPosPtr;
129 : Double *__restrict__ convFunc=(*(convFuncStore_p.rdata)).getStorage(Dummy);
130 : Double * __restrict__ freq=vbs.freq_p.getStorage(Dummy);
131 : Bool * __restrict__ rowFlag=vbs.rowFlag_p.getStorage(Dummy);
132 :
133 : Float * __restrict__ imagingWeight = vbs.imagingWeight_p.getStorage(Dummy);
134 : Double * __restrict__ uvw = vbs.uvw_p.getStorage(Dummy);
135 : Bool * __restrict__ flagCube = vbs.flagCube_p.getStorage(Dummy);
136 : Complex * __restrict__ visCube = vbs.visCube_p.getStorage(Dummy);
137 : Double * __restrict__ scale = uvwScale_p.getStorage(Dummy);
138 : Double * __restrict__ offset = offset_p.getStorage(Dummy);
139 : Float * __restrict__ samplingPtr = sampling.getStorage(Dummy);
140 : Double * __restrict__ posPtr=pos.getStorage(Dummy);
141 : Int * __restrict__ locPtr=loc.getStorage(Dummy);
142 : Int * __restrict__ offPtr=off.getStorage(Dummy);
143 : Double * __restrict__ sumwtPtr = sumwt.getStorage(Dummy);
144 : Int * __restrict__ ilocPtr=iloc.getStorage(Dummy);
145 : Int * __restrict__ supportPtr = support.getStorage(Dummy);
146 : Int nDim = vbs.uvw_p.shape()[0];
147 :
148 : // cacheAxisIncrements(nx,ny,nGridPol, nGridChan);
149 : cacheAxisIncrements(grid.shape().asVector());
150 :
151 : for(Int irow=rbeg; irow < rend; irow++){ // For all rows
152 :
153 : if(!rowFlag[irow]){ // If the row is not flagged
154 :
155 : for(Int ichan=0; ichan< nDataChan; ichan++){ // For all channels
156 :
157 : // if (vbs.imagingWeight(ichan,irow)!=0.0) { // If weights are not zero
158 : if (imagingWeight[ichan+irow*nDataChan]!=0.0) { // If weights are not zero
159 : achan=chanMap_p(ichan);
160 :
161 : if((achan>=0) && (achan<nGridChan)) { // If selected channels are valid
162 :
163 : // sgrid(pos,loc,off, phasor, irow,
164 : // vbs.uvw,dphase_p[irow], vbs.freq[ichan],
165 : // uvwScale_p, offset_p, sampling);
166 : sgrid(nDim,posPtr,locPtr,offPtr, phasor, irow,
167 : uvw,dphase_p[irow], freq[ichan],
168 : scale, offset, samplingPtr);
169 :
170 : if (onGrid(nx, ny, loc, support)) { // If the data co-ords. are with-in the grid
171 :
172 : for(Int ipol=0; ipol< nDataPol; ipol++) { // For all polarizations
173 : // if((!vbs.flagCube(ipol,ichan,irow))){ // If the pol. & chan. specific
174 : if((!flagCube[ipol+ichan*nDataPol+irow*nDataChan*nDataPol])){
175 : apol=polMap_p(ipol);
176 : if ((apol>=0) && (apol<nGridPol)) {
177 : // igrdpos(2)=apol; igrdpos(3)=achan;
178 : iPosPtr[2]=apol; iPosPtr[3]=achan;
179 : norm=0.0;
180 :
181 : imgWt=imagingWeight[ichan+irow*nDataChan];
182 :
183 : Int idopsf = (dopsf==true);
184 : Int ncf=(*convFuncStore_p.rdata).shape()(0);
185 : Int nrow = vbs.nRow();
186 : gridengine(&nx, &ny, &nGridPol, &nGridChan, (double*)gridStore,
187 : &ncf, (double*)convFunc, samplingPtr, supportPtr,
188 : locPtr, offPtr, &achan, &apol, &norm, (float *)&phasor,
189 : &ipol, &ichan, &irow, imagingWeight, &idopsf,
190 : (float *)visCube, &nDataPol, &nDataChan,&nrow);
191 :
192 : // if(dopsf) nvalue=Complex(imgWt);
193 : // else nvalue=imgWt*
194 : // // (vbs.visCube(ipol,ichan,irow)*phasor);
195 : // (visCube[ipol+ichan*nDataPol+irow*nDataPol*nDataChan]*phasor);
196 :
197 : // for(Int iy=-supportPtr[1]; iy <= supportPtr[1]; iy++)
198 : // {
199 : // ilocPtr[1]=abs((int)(samplingPtr[1]*iy+offPtr[1]));
200 : // // igrdpos(1)=loc(1)+iy;
201 : // iPosPtr[1]=locPtr[1]+iy;
202 : // // wt = convFunc[ilocPtr[1]];
203 : // for(Int ix=-supportPtr[0]; ix <= supportPtr[0]; ix++)
204 : // {
205 : // ilocPtr[0]=abs((int)(samplingPtr[0]*ix+offPtr[0]));
206 : // wt = convFunc[iloc[0]]*convFunc[iloc[1]];
207 : // // wt *= convFunc[ilocPtr[0]];
208 :
209 : // //igrdpos(0)=loc(0)+ix;
210 : // iPosPtr[0]=locPtr[0]+ix;
211 :
212 :
213 : // // grid(grdpos) += nvalue*wt;
214 :
215 : // // gridStore[iPosPtr[0] +
216 : // // iPosPtr[1]*incPtr_p[1] +
217 : // // iPosPtr[2]*incPtr_p[2] +
218 : // // iPosPtr[3]*incPtr_p[3]].real() += (nvalue.real()*wt);
219 : // // gridStore[iPosPtr[0] +
220 : // // iPosPtr[1]*incPtr_p[1] +
221 : // // iPosPtr[2]*incPtr_p[2] +
222 : // // iPosPtr[3]*incPtr_p[3]].imag() += (nvalue.imag()*wt);
223 :
224 : // // The following uses raw index on the 4D grid
225 : // // addTo4DArray(gridStore,iPosPtr,nvalue,wt);
226 :
227 : // addtogrid((double *)gridStore, iPosPtr, (float *)&nvalue, (double *)&wt,
228 : // &nx, &ny, &nGridPol, &nGridChan);
229 :
230 : // norm+=wt;
231 : // }
232 : // }
233 : sumwtPtr[apol+achan*nGridPol]+=imgWt*norm;
234 : // sumwt(apol,achan)+=imgWt*norm;
235 : }
236 : }
237 : }
238 : }
239 : }
240 : }
241 : }
242 : }
243 : }
244 : // T *tt=(T *)gridStore;
245 : // grid.putStorage(tt,gDummy);
246 : }
247 : */
248 : //
249 : //-----------------------------------------------------------------------------------
250 : // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding)
251 : //
252 0 : void VisibilityResampler::GridToData(VBStore& vbs, const Array<Complex>& grid)
253 : // void VisibilityResampler::GridToData(VBStore& vbs, Array<Complex>& grid)
254 : {
255 : using casacore::operator*;
256 :
257 : Int nDataChan, nDataPol, nGridPol, nGridChan, nx, ny;
258 : Int achan, apol, rbeg, rend;
259 0 : Vector<Float> sampling(2);
260 0 : Vector<Int> support(2),loc(2), off(2), iloc(2);
261 0 : Vector<Double> pos(2);
262 :
263 0 : IPosition grdpos(4);
264 :
265 0 : Double norm=0, wt;
266 0 : Complex phasor, nvalue;
267 :
268 : // rbeg=0;
269 : // rend=vbs.nRow_p;
270 0 : rbeg = vbs.beginRow_p;
271 0 : rend = vbs.endRow_p;
272 : // cerr << rbeg << " " << rend << " " << vbs.nRow() << endl;
273 0 : nx = grid.shape()[0]; ny = grid.shape()[1];
274 0 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
275 :
276 0 : nDataPol = vbs.flagCube_p.shape()[0];
277 0 : nDataChan = vbs.flagCube_p.shape()[1];
278 :
279 0 : sampling[0] = sampling[1] = convFuncStore_p.sampling[0];
280 0 : support(0) = convFuncStore_p.xSupport[0];
281 0 : support(1) = convFuncStore_p.ySupport[0];
282 :
283 : Bool Dummy,vbcDummy;
284 0 : const Complex *__restrict__ gridStore = grid.getStorage(Dummy);
285 0 : Vector<Int> igrdpos(4);
286 0 : Int *__restrict__ iPosPtr = igrdpos.getStorage(Dummy);
287 0 : const Int *__restrict__ iPosPtrConst = iPosPtr;
288 0 : Double *__restrict__ convFunc=(*(convFuncStore_p.rdata)).getStorage(Dummy);
289 0 : Double *__restrict__ freq=vbs.freq_p.getStorage(Dummy);
290 0 : Bool *__restrict__ rowFlag=vbs.rowFlag_p.getStorage(Dummy);
291 :
292 : //UNUSED: Float * __restrict__ imagingWeight = vbs.imagingWeight_p.getStorage(Dummy);
293 0 : Double * __restrict__ uvw = vbs.uvw_p.getStorage(Dummy);
294 0 : Bool * __restrict__ flagCube = vbs.flagCube_p.getStorage(Dummy);
295 0 : Complex * __restrict__ visCube = vbs.visCube_p.getStorage(vbcDummy);
296 0 : Double * __restrict__ scale = uvwScale_p.getStorage(Dummy);
297 0 : Double * __restrict__ offset = offset_p.getStorage(Dummy);
298 0 : Int * __restrict__ supportPtr = support.getStorage(Dummy);
299 0 : Float * __restrict__ samplingPtr = sampling.getStorage(Dummy);
300 0 : Double * __restrict__ posPtr=pos.getStorage(Dummy);
301 0 : Int * __restrict__ locPtr=loc.getStorage(Dummy);
302 0 : Int * __restrict__ offPtr=off.getStorage(Dummy);
303 0 : Int * __restrict__ ilocPtr=iloc.getStorage(Dummy);
304 0 : Int nDim = vbs.uvw_p.shape()(0);
305 :
306 : // cacheAxisIncrements(nx,ny,nGridPol, nGridChan);
307 0 : cacheAxisIncrements(grid.shape().asVector());
308 : //UNUSED: Int xSupport_l, ySupport_l;
309 : //UNUSED: Float sampling_l;
310 0 : for(Int irow=rbeg; irow < rend; irow++) {
311 0 : if(!rowFlag[irow]) {
312 0 : for (Int ichan=0; ichan < nDataChan; ichan++) {
313 0 : achan=chanMap_p(ichan);
314 :
315 0 : if((achan>=0) && (achan<nGridChan)) {
316 : // sgrid(pos,loc,off,phasor,irow,vbs.uvw,
317 : // dphase_p[irow],vbs.freq[ichan],
318 : // uvwScale_p,offset_p,sampling);
319 :
320 0 : sgrid(nDim,posPtr,locPtr,offPtr,phasor,irow,uvw,
321 0 : dphase_p[irow],freq[ichan],
322 : scale,offset,samplingPtr);
323 :
324 0 : if (onGrid(nx, ny, loc, support)) {
325 0 : for(Int ipol=0; ipol < nDataPol; ipol++) {
326 :
327 0 : if(!flagCube[ipol+ichan*nDataPol+irow*nDataChan*nDataPol]) {
328 0 : apol=polMap_p(ipol);
329 :
330 0 : if((apol>=0) && (apol<nGridPol)) {
331 : // igrdpos(2)=apol; igrdpos(3)=achan;
332 0 : iPosPtr[2]=apol; iPosPtr[3]=achan;
333 0 : nvalue=0.0;
334 0 : norm=0.0;
335 :
336 0 : for(Int iy=-supportPtr[1]; iy <= supportPtr[1]; iy++)
337 : {
338 0 : ilocPtr[1]=abs((Int)(samplingPtr[1]*iy+offPtr[1]));
339 : // igrdpos(1)=loc(1)+iy;
340 0 : iPosPtr[1]=locPtr[1]+iy;
341 : // wt = convFunc[ilocPtr[1]];
342 0 : for(Int ix=-supportPtr[0]; ix <= supportPtr[0]; ix++)
343 : {
344 0 : ilocPtr[0]=abs((Int)(samplingPtr[0]*ix+offPtr[0]));
345 : // igrdpos(0)=loc(0)+ix;
346 0 : iPosPtr[0]=locPtr[0]+ix;
347 0 : wt=convFunc[ilocPtr[0]]*convFunc[ilocPtr[1]];
348 : // wt *= convFunc[ilocPtr[0]];
349 0 : norm+=wt;
350 : // nvalue+=wt*grid(grdpos);
351 : // The following uses raw index on the 4D grid
352 0 : nvalue+=wt*getFrom4DArray(gridStore,iPosPtrConst);
353 : }
354 : }
355 0 : visCube[ipol+ichan*nDataPol+irow*nDataChan*nDataPol]=(nvalue*conj(phasor))/norm;
356 : }
357 : }
358 : }
359 : }
360 : }
361 : }
362 : }
363 : }
364 0 : Complex *tt=(Complex *) visCube;
365 0 : vbs.visCube_p.putStorage(tt,vbcDummy);
366 0 : }
367 : //
368 : //-----------------------------------------------------------------------------------
369 : //
370 : // void VisibilityResampler::sgrid(Int& uvwDim,Double* __restrict__ pos,
371 : // Int* __restrict__ loc,
372 : // Int* __restrict__ off,
373 : // Complex& phasor, const Int& irow,
374 : // // const Matrix<Double>& __restrict__ uvw,
375 : // const Double* __restrict__ uvw,
376 : // const Double& __restrict__ dphase,
377 : // const Double& __restrict__ freq,
378 : // const Double* __restrict__ scale,
379 : // const Double* __restrict__ offset,
380 : // const Float* __restrict__ sampling) __restrict__
381 : // // const Vector<Double>& __restrict__ scale,
382 : // // const Vector<Double>& __restrict__ offset,
383 : // // const Vector<Float>& __restrict__ sampling) __restrict__
384 : // {
385 : // Double phase;
386 : // // Int ndim=pos.shape()(0);
387 : // Int ndim=2;
388 :
389 : // for(Int idim=0;idim<ndim;idim++)
390 : // {
391 : // pos[idim]=scale[idim]*uvw[idim+irow*uvwDim]*freq/C::c+offset[idim];
392 : // loc[idim]=(Int)std::floor(pos[idim]+0.5);
393 : // off[idim]=(Int)std::floor(((loc[idim]-pos[idim])*sampling[idim])+0.5);
394 : // }
395 :
396 : // if (dphase != 0.0)
397 : // {
398 : // phase=-2.0*C::pi*dphase*freq/C::c;
399 : // phasor=Complex(cos(phase), sin(phase));
400 : // }
401 : // else
402 : // phasor=Complex(1.0);
403 : // }
404 0 : void VisibilityResampler::ComputeResiduals(VBStore& vbs)
405 : {
406 0 : Int rbeg = vbs.beginRow_p, rend = vbs.endRow_p;
407 0 : IPosition vbDataShape=vbs.modelCube_p.shape();
408 0 : IPosition start(vbDataShape), last(vbDataShape);
409 0 : start=0; start(2)=rbeg;
410 0 : last(2)=rend; //last=last-1;
411 :
412 0 : if (!vbs.useCorrected_p)
413 : {
414 0 : for(uInt ichan = start(0); ichan < last(0); ichan++)
415 0 : for(uInt ipol = start(1); ipol < last(1); ipol++)
416 0 : for(uInt irow = start(2); irow < last(2); irow++)
417 0 : vbs.modelCube_p(ichan,ipol,irow) = vbs.modelCube_p(ichan,ipol,irow) - vbs.visCube_p(ichan,ipol,irow);
418 : // vbs.modelCube_p(ichan,ipol,irow) = vbs.visCube_p(ichan,ipol,irow) - vbs.modelCube_p(ichan,ipol,irow);
419 : }
420 : else
421 : {
422 0 : for(uInt ichan = start(0); ichan < last(0); ichan++)
423 0 : for(uInt ipol = start(1); ipol < last(1); ipol++)
424 0 : for(uInt irow = start(2); irow < last(2); irow++)
425 0 : vbs.modelCube_p(ichan,ipol,irow) = vbs.modelCube_p(ichan,ipol,irow) - vbs.correctedCube_p(ichan,ipol,irow);
426 : //vbs.modelCube_p(ichan,ipol,irow) = vbs.correctedCube_p(ichan,ipol,irow) - vbs.modelCube_p(ichan,ipol,irow);
427 : }
428 :
429 : // Slicer mySlice(start,last,Slicer::endIsLast);
430 : // Cube<Complex> slicedModel, slicedData, slicedCorrected;
431 : // if (!vbs.useCorrected_p)
432 : // {
433 : // {
434 : // async::MutexLocker tt(*myMutex_p);
435 : // slicedModel = Cube<Complex>(vbs.modelCube_p(mySlice));
436 : // slicedData = Cube<Complex>(vbs.visCube_p(mySlice));
437 : // }
438 : // slicedModel -= slicedData;
439 : // }
440 : // else
441 : // {
442 : // {
443 : // async::MutexLocker tt(*myMutex_p);
444 : // slicedModel = Cube<Complex>(vbs.modelCube_p(mySlice));
445 : // slicedCorrected = Cube<Complex>(vbs.correctedCube_p(mySlice));
446 : // }
447 : // slicedModel -= slicedCorrected;
448 : // }
449 0 : }
450 :
451 :
452 :
453 :
454 : #define DoOld 1
455 :
456 : template <class T>
457 0 : void VisibilityResampler::DataToGridImpl_p(Array<T>& grid, VBStore& vbs, const Bool& dopsf,
458 : Matrix<Double>& sumwt,
459 : Bool /*useConjFreqCF*/)
460 : {
461 : using casacore::operator*;
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 : 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 : using namespace casacore;
695 : };// end namespace casa
|
