Line data Source code
1 : //# WProjectFT.cc: Implementation of WProjectFT class
2 : //# Copyright (C) 2003-2016
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <msvis/MSVis/VisibilityIterator2.h>
29 : #include <casacore/casa/Quanta/UnitMap.h>
30 : #include <casacore/casa/Quanta/MVTime.h>
31 : #include <casacore/casa/Quanta/UnitVal.h>
32 : #include <casacore/casa/Utilities/CountedPtr.h>
33 : #include <casacore/measures/Measures/Stokes.h>
34 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
35 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
36 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
37 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
38 : #include <casacore/coordinates/Coordinates/Projection.h>
39 : #include <casacore/ms/MeasurementSets/MSColumns.h>
40 : #include <casacore/casa/BasicSL/Constants.h>
41 : #include <casacore/scimath/Mathematics/FFTServer.h>
42 : #include <synthesis/TransformMachines2/WProjectFT.h>
43 : #include <synthesis/TransformMachines2/WPConvFunc.h>
44 : #include <casacore/scimath/Mathematics/RigidVector.h>
45 : #include <msvis/MSVis/StokesVector.h>
46 : #include <synthesis/TransformMachines/StokesImageUtil.h>
47 : #include <casacore/images/Images/ImageInterface.h>
48 : #include <casacore/images/Images/PagedImage.h>
49 : #include <casacore/casa/Containers/Block.h>
50 : #include <casacore/casa/Containers/Record.h>
51 : #include <casacore/casa/Arrays/ArrayLogical.h>
52 : #include <casacore/casa/Arrays/ArrayMath.h>
53 : #include <casacore/casa/Arrays/Array.h>
54 : #include <casacore/casa/Arrays/MaskedArray.h>
55 : #include <casacore/casa/Arrays/Vector.h>
56 : #include <casacore/casa/Arrays/Slice.h>
57 : #include <casacore/casa/Arrays/Matrix.h>
58 : #include <casacore/casa/Arrays/Cube.h>
59 : #include <casacore/casa/Arrays/MatrixIter.h>
60 : #include <casacore/casa/BasicSL/String.h>
61 : #include <casacore/casa/Utilities/Assert.h>
62 : #include <casacore/casa/Exceptions/Error.h>
63 : #include <iostream>
64 : #include <iomanip>
65 : #include <casacore/lattices/Lattices/ArrayLattice.h>
66 : #include <casacore/lattices/Lattices/SubLattice.h>
67 : #include <casacore/lattices/LRegions/LCBox.h>
68 : #include <casacore/lattices/LEL/LatticeExpr.h>
69 : #include <casacore/lattices/Lattices/LatticeCache.h>
70 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
71 : #include <casacore/lattices/Lattices/LatticeIterator.h>
72 : #include <casacore/lattices/Lattices/LatticeStepper.h>
73 : #include <casacore/casa/Utilities/CompositeNumber.h>
74 : #include <casacore/casa/OS/Timer.h>
75 : #include <casacore/casa/OS/HostInfo.h>
76 : #include <sstream>
77 : #ifdef _OPENMP
78 : #include <omp.h>
79 : #endif
80 : using namespace casacore;
81 : namespace casa { //# NAMESPACE CASA - BEGIN
82 : namespace refim { //#namespace for imaging refactoring
83 :
84 : using namespace casacore;
85 : using namespace casa;
86 : using namespace casacore;
87 : using namespace casa::vi;
88 : using namespace casacore;
89 : using namespace casa::refim;
90 :
91 0 : WProjectFT::WProjectFT( Int nWPlanes, Long icachesize, Int itilesize,
92 0 : Bool usezero, Bool useDoublePrec, const Double minW, const Double maxW, const Double rmsW)
93 : : FTMachine(), padding_p(1.0), nWPlanes_p(nWPlanes),
94 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
95 : gridder(0), isTiled(false),
96 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)), usezero_p(usezero),
97 0 : machineName_p("WProjectFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), minW_p(minW), maxW_p(maxW), rmsW_p(rmsW)
98 : {
99 0 : convSize=0;
100 0 : tangentSpecified_p=false;
101 0 : lastIndex_p=0;
102 0 : useDoubleGrid_p=useDoublePrec;
103 0 : wpConvFunc_p=new WPConvFunc(minW, maxW, rmsW);
104 :
105 0 : }
106 78 : WProjectFT::WProjectFT(Int nWPlanes,
107 : MPosition mLocation,
108 : Long icachesize, Int itilesize,
109 78 : Bool usezero, Float padding, Bool useDoublePrec, const Double minW, const Double maxW, const Double rmsW)
110 : : FTMachine(), padding_p(padding), nWPlanes_p(nWPlanes),
111 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
112 : gridder(0), isTiled(false),
113 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
114 : usezero_p(usezero),
115 78 : machineName_p("WProjectFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), minW_p(minW), maxW_p(maxW), rmsW_p(rmsW)
116 :
117 : {
118 78 : convSize=0;
119 78 : savedWScale_p=0.0;
120 78 : tangentSpecified_p=false;
121 78 : mLocation_p=mLocation;
122 78 : lastIndex_p=0;
123 78 : wpConvFunc_p=new WPConvFunc(minW, maxW, rmsW);
124 78 : useDoubleGrid_p=useDoublePrec;
125 78 : }
126 2 : WProjectFT::WProjectFT(
127 : Int nWPlanes, MDirection mTangent,
128 : MPosition mLocation,
129 : Long icachesize, Int itilesize,
130 2 : Bool usezero, Float padding, Bool useDoublePrec, const Double minW, const Double maxW, const Double rmsW)
131 : : FTMachine(), padding_p(padding), nWPlanes_p(nWPlanes),
132 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
133 : gridder(0), isTiled(false),
134 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
135 : usezero_p(usezero),
136 2 : machineName_p("WProjectFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), minW_p(minW), maxW_p(maxW), rmsW_p(rmsW)
137 : {
138 2 : convSize=0;
139 2 : savedWScale_p=0.0;
140 2 : mTangent_p=mTangent;
141 2 : tangentSpecified_p=true;
142 2 : mLocation_p=mLocation;
143 2 : lastIndex_p=0;
144 2 : wpConvFunc_p=new WPConvFunc(minW, maxW, rmsW);
145 2 : useDoubleGrid_p=useDoublePrec;
146 2 : }
147 :
148 0 : WProjectFT::WProjectFT(const RecordInterface& stateRec)
149 0 : : FTMachine(),machineName_p("WProjectFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0)
150 : {
151 : // Construct from the input state record
152 0 : String error;
153 0 : if (!fromRecord(error, stateRec)) {
154 0 : throw (AipsError("Failed to create WProjectFT: " + error));
155 : };
156 0 : }
157 : //----------------------------------------------------------------------
158 30 : WProjectFT& WProjectFT::operator=(const WProjectFT& other)
159 : {
160 30 : if(this!=&other) {
161 : //Do the base parameters
162 30 : FTMachine::operator=(other);
163 :
164 :
165 30 : padding_p=other.padding_p;
166 30 : nWPlanes_p=other.nWPlanes_p;
167 30 : imageCache=other.imageCache;
168 30 : cachesize=other.cachesize;
169 30 : tilesize=other.tilesize;
170 30 : if(other.gridder==0)
171 30 : gridder=0;
172 : else{
173 0 : uvScale.resize();
174 0 : uvOffset.resize();
175 0 : uvScale=other.uvScale;
176 0 : uvOffset=other.uvOffset;
177 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
178 0 : uvScale, uvOffset,
179 0 : "SF");
180 : }
181 :
182 30 : isTiled=other.isTiled;
183 : //lattice=other.lattice;
184 : //arrayLattice=other.arrayLattice;
185 30 : lattice=0;
186 30 : arrayLattice=0;
187 :
188 30 : maxAbsData=other.maxAbsData;
189 30 : centerLoc=other.centerLoc;
190 30 : offsetLoc=other.offsetLoc;
191 30 : usezero_p=other.usezero_p;
192 30 : machineName_p=other.machineName_p;
193 30 : wpConvFunc_p=other.wpConvFunc_p;
194 30 : timemass_p=0.0;
195 30 : timegrid_p=0.0;
196 30 : timedegrid_p=0.0;
197 30 : minW_p=other.minW_p;
198 30 : maxW_p=other.maxW_p;
199 30 : rmsW_p=other.rmsW_p;
200 : };
201 30 : return *this;
202 : };
203 :
204 : //----------------------------------------------------------------------
205 30 : WProjectFT::WProjectFT(const WProjectFT& other) :FTMachine(), machineName_p("WProjectFT"),timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0)
206 : {
207 30 : operator=(other);
208 30 : }
209 :
210 30 : FTMachine* WProjectFT::cloneFTM(){
211 30 : return new WProjectFT(*this);
212 : }
213 :
214 : //----------------------------------------------------------------------
215 86 : void WProjectFT::init() {
216 : /* if((padding_p*padding_p*image->shape().product())>cachesize) {
217 : isTiled=true;
218 : nx = image->shape()(0);
219 : ny = image->shape()(1);
220 : npol = image->shape()(2);
221 : nchan = image->shape()(3);
222 : }
223 : else {*/
224 : // We are padding.
225 86 : isTiled=false;
226 172 : CompositeNumber cn(uInt(image->shape()(0)*2));
227 86 : nx = cn.nextLargerEven(Int(padding_p*Float(image->shape()(0))-0.5));
228 86 : ny = cn.nextLargerEven(Int(padding_p*Float(image->shape()(1))-0.5));
229 86 : npol = image->shape()(2);
230 86 : nchan = image->shape()(3);
231 : //}
232 :
233 : // if(image->shape().product()>cachesize) {
234 : // isTiled=true;
235 : // }
236 : // else {
237 : // isTiled=false;
238 : // }
239 : //The Tiled version need some fixing: sof or now
240 86 : isTiled=false;
241 :
242 :
243 86 : sumWeight.resize(npol, nchan);
244 :
245 86 : wConvSize=max(0, nWPlanes_p);
246 86 : convSupport.resize(wConvSize);
247 86 : convSupport=0;
248 :
249 86 : uvScale.resize(3);
250 86 : uvScale=0.0;
251 86 : uvScale(0)=Float(nx)*image->coordinates().increment()(0);
252 86 : uvScale(1)=Float(ny)*image->coordinates().increment()(1);
253 86 : if(savedWScale_p==0.0){
254 50 : uvScale(2)=Float(wConvSize)*abs(image->coordinates().increment()(0));
255 : }
256 : else{
257 36 : uvScale(2)=savedWScale_p;
258 : }
259 86 : uvOffset.resize(3);
260 86 : uvOffset(0)=nx/2;
261 86 : uvOffset(1)=ny/2;
262 86 : uvOffset(2)=0;
263 :
264 :
265 :
266 86 : if(gridder) delete gridder; gridder=0;
267 172 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
268 86 : uvScale, uvOffset,
269 86 : "SF");
270 :
271 : // Set up image cache needed for gridding.
272 86 : if(imageCache) delete imageCache; imageCache=0;
273 :
274 : // The tile size should be large enough that the
275 : // extended convolution function can fit easily
276 86 : if(isTiled) {
277 0 : Float tileOverlap=0.5;
278 0 : tilesize=min(256,tilesize);
279 0 : IPosition tileShape=IPosition(4,tilesize,tilesize,npol,nchan);
280 0 : Vector<Float> tileOverlapVec(4);
281 0 : tileOverlapVec=0.0;
282 0 : tileOverlapVec(0)=tileOverlap;
283 0 : tileOverlapVec(1)=tileOverlap;
284 0 : Int tmpCacheVal=static_cast<Int>(cachesize);
285 0 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
286 : tileOverlapVec,
287 0 : (tileOverlap>0.0));
288 :
289 : }
290 86 : }
291 :
292 : // This is nasty, we should use CountedPointers here.
293 220 : WProjectFT::~WProjectFT() {
294 110 : if(imageCache) delete imageCache; imageCache=0;
295 110 : if(gridder) delete gridder; gridder=0;
296 : /*
297 : if(arrayLattice) delete arrayLattice; arrayLattice=0;
298 :
299 : Int numofmodels=convFunctions_p.nelements();
300 : for (Int k=0; k< numofmodels; ++k){
301 : delete convFunctions_p[k];
302 : delete convSupportBlock_p[k];
303 :
304 : }
305 : */
306 : // convFuctions_p.resize();
307 : // convSupportBlock_p.resize();
308 :
309 220 : }
310 :
311 :
312 40 : void WProjectFT::setConvFunc(CountedPtr<WPConvFunc>& pbconvFunc){
313 :
314 40 : wpConvFunc_p=pbconvFunc;
315 40 : }
316 40 : CountedPtr<WPConvFunc>& WProjectFT::getConvFunc(){
317 :
318 40 : return wpConvFunc_p;
319 : }
320 :
321 1544 : void WProjectFT::findConvFunction(const ImageInterface<Complex>& image,
322 : const VisBuffer2& vb) {
323 :
324 :
325 :
326 :
327 :
328 1544 : wpConvFunc_p->findConvFunction(image, vb, wConvSize, uvScale, uvOffset,
329 1544 : padding_p,
330 1544 : convSampling,
331 1544 : convFunc, convSize, convSupport,
332 1544 : savedWScale_p);
333 :
334 1544 : nWPlanes_p=convSupport.nelements();
335 1544 : wConvSize=max(1,nWPlanes_p);
336 1544 : uvScale(2)=savedWScale_p;
337 :
338 :
339 :
340 1544 : }
341 :
342 13 : void WProjectFT::initializeToVis(ImageInterface<Complex>& iimage,
343 : const VisBuffer2& vb)
344 : {
345 13 : image=&iimage;
346 13 : toVis_p=true;
347 13 : ok();
348 :
349 : // if(convSize==0) {
350 13 : init();
351 : // }
352 13 : findConvFunction(*image, vb);
353 :
354 :
355 : // Initialize the maps for polarization and channel. These maps
356 : // translate visibility indices into image indices
357 13 : initMaps(vb);
358 :
359 : // nx = image->shape()(0);
360 : // ny = image->shape()(1);
361 : // npol = image->shape()(2);
362 : // nchan = image->shape()(3);
363 :
364 :
365 13 : isTiled=false;
366 : // If we are memory-based then read the image in and create an
367 : // ArrayLattice otherwise just use the PagedImage
368 : /*if(isTiled) {
369 : lattice=CountedPtr<Lattice<Complex> > (image, false);
370 : }
371 : else {
372 : }
373 : */
374 : //AlwaysAssert(lattice, AipsError);
375 :
376 13 : prepGridForDegrid();
377 13 : }
378 :
379 13 : void WProjectFT::prepGridForDegrid(){
380 26 : IPosition gridShape(4, nx, ny, npol, nchan);
381 13 : griddedData.resize(gridShape);
382 13 : griddedData=Complex(0.0);
383 :
384 :
385 26 : IPosition stride(4, 1);
386 26 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
387 39 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
388 39 : IPosition trc(blc+image->shape()-stride);
389 :
390 26 : IPosition start(4, 0);
391 13 : griddedData(blc, trc) = image->getSlice(start, image->shape());
392 :
393 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
394 13 : arrayLattice = new ArrayLattice<Complex>(griddedData);
395 13 : lattice=arrayLattice;
396 :
397 13 : logIO() << LogIO::DEBUGGING << "Starting FFT of image" << LogIO::POST;
398 :
399 :
400 26 : Vector<Float> sincConvX(nx);
401 16173 : for (Int ix=0;ix<nx;ix++) {
402 16160 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
403 16160 : if(ix==nx/2) {
404 13 : sincConvX(ix)=1.0;
405 : }
406 : else {
407 16147 : sincConvX(ix)=sin(x)/x;
408 : }
409 : }
410 26 : Vector<Float> sincConvY(ny);
411 16173 : for (Int ix=0;ix<ny;ix++) {
412 16160 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
413 16160 : if(ix==ny/2) {
414 13 : sincConvY(ix)=1.0;
415 : }
416 : else {
417 16147 : sincConvY(ix)=sin(x)/x;
418 : }
419 : }
420 :
421 26 : Vector<Complex> correction(nx);
422 13 : correction=Complex(1.0, 0.0);
423 : // Do the Grid-correction
424 26 : IPosition cursorShape(4, nx, 1, 1, 1);
425 26 : IPosition axisPath(4, 0, 1, 2, 3);
426 26 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
427 26 : LatticeIterator<Complex> lix(*lattice, lsx);
428 22453 : for(lix.reset();!lix.atEnd();lix++) {
429 22440 : Int iy=lix.position()(1);
430 22440 : gridder->correctX1D(correction,iy);
431 34250640 : for (Int ix=0;ix<nx;ix++) {
432 34228200 : correction(ix)/=(sincConvX(ix)*sincConvY(iy));
433 : }
434 22440 : lix.rwVectorCursor()/=correction;
435 : }
436 :
437 : // Now do the FFT2D in place
438 13 : LatticeFFT::cfft2d(*lattice);
439 :
440 13 : logIO() << LogIO::DEBUGGING << "Finished FFT" << LogIO::POST;
441 :
442 :
443 :
444 13 : }
445 :
446 :
447 13 : void WProjectFT::finalizeToVis()
448 : {
449 13 : logIO() << LogOrigin("WProjectFT", "finalizeToVis") << LogIO::NORMAL;
450 13 : logIO() <<LogIO::NORMAL2<< "Time to degrid " << timedegrid_p << LogIO::POST ;
451 13 : timedegrid_p=0.0;
452 13 : if(!arrayLattice.null()) arrayLattice=0;
453 13 : if(!lattice.null()) lattice=0;
454 13 : griddedData.resize();
455 13 : if(isTiled) {
456 :
457 :
458 :
459 0 : AlwaysAssert(imageCache, AipsError);
460 0 : AlwaysAssert(image, AipsError);
461 0 : ostringstream o;
462 0 : imageCache->flush();
463 0 : imageCache->showCacheStatistics(o);
464 0 : logIO() << o.str() << LogIO::POST;
465 : }
466 13 : }
467 :
468 :
469 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
470 : // grid.
471 73 : void WProjectFT::initializeToSky(ImageInterface<Complex>& iimage,
472 : Matrix<Float>& weight,
473 : const VisBuffer2& vb)
474 : {
475 : // image always points to the image
476 73 : image=&iimage;
477 73 : toVis_p=false;
478 :
479 : // if(convSize==0) {
480 73 : init();
481 : // }
482 73 : findConvFunction(*image, vb);
483 :
484 :
485 : // Initialize the maps for polarization and channel. These maps
486 : // translate visibility indices into image indices
487 73 : initMaps(vb);
488 :
489 : // nx = image->shape()(0);
490 : // ny = image->shape()(1);
491 : // npol = image->shape()(2);
492 : // nchan = image->shape()(3);
493 :
494 : // if(image->shape().product()>cachesize) {
495 : // isTiled=true;
496 : // }
497 : // else {
498 : // isTiled=false;
499 : // }
500 73 : isTiled=false;
501 73 : sumWeight=0.0;
502 73 : weight.resize(sumWeight.shape());
503 73 : weight=0.0;
504 :
505 : // Initialize for in memory or to disk gridding. lattice will
506 : // point to the appropriate Lattice, either the ArrayLattice for
507 : // in memory gridding or to the image for to disk gridding.
508 73 : if(isTiled) {
509 0 : imageCache->flush();
510 0 : image->set(Complex(0.0));
511 0 : lattice=CountedPtr<Lattice<Complex> > (image, false);
512 : }
513 : else {
514 146 : IPosition gridShape(4, nx, ny, npol, nchan);
515 73 : if(!useDoubleGrid_p){
516 0 : griddedData.resize(gridShape);
517 0 : griddedData=Complex(0.0);
518 : }
519 : else{
520 : //griddedData.resize();
521 73 : griddedData2.resize(gridShape);
522 73 : griddedData2=DComplex(0.0);
523 : }
524 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
525 :
526 : }
527 : //AlwaysAssert(lattice, AipsError);
528 :
529 73 : }
530 :
531 73 : void WProjectFT::finalizeToSky()
532 : {
533 73 : logIO() << LogOrigin("WProjectFT", "finalizeToSky") << LogIO::NORMAL;
534 73 : logIO() <<LogIO::NORMAL2<< "Time to massage data " << timemass_p << LogIO::POST;
535 73 : logIO() << LogIO::NORMAL2 <<"Time to grid data " << timegrid_p << LogIO::POST;
536 73 : timemass_p=0.0;
537 73 : timegrid_p=0.0;
538 : // Now we flush the cache and report statistics
539 : // For memory based, we don't write anything out yet.
540 73 : if(isTiled) {
541 0 : logIO() << LogOrigin("WProjectFT", "finalizeToSky") << LogIO::NORMAL;
542 :
543 0 : AlwaysAssert(image, AipsError);
544 0 : AlwaysAssert(imageCache, AipsError);
545 0 : imageCache->flush();
546 0 : ostringstream o;
547 0 : imageCache->showCacheStatistics(o);
548 0 : logIO() << o.str() << LogIO::POST;
549 : }
550 73 : }
551 :
552 0 : Array<Complex>* WProjectFT::getDataPointer(const IPosition& centerLoc2D,
553 : Bool readonly) {
554 : Array<Complex>* result;
555 : // Is tiled: get tiles and set up offsets
556 0 : centerLoc(0)=centerLoc2D(0);
557 0 : centerLoc(1)=centerLoc2D(1);
558 0 : result=&imageCache->tile(offsetLoc, centerLoc, readonly);
559 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
560 0 : return result;
561 : }
562 :
563 : #define NEED_UNDERSCORES
564 : #if defined(NEED_UNDERSCORES)
565 : #define gwgrid gwgrid_
566 : #define gwproj gwproj_
567 : #define dwproj dwproj_
568 : #define sectgwgridd sectgwgridd_
569 : #define sectgwgrids sectgwgrids_
570 : #define sectdwgrid sectdwgrid_
571 : #define locuvw locuvw_
572 : #endif
573 :
574 : extern "C" {
575 : void locuvw(const Double*, const Double*, const Double*, const Int*, const Double*, const Double*, const Int*,
576 : Int*, Int*, Complex*, const Int*, const Int*, const Double*);
577 : //Double precision gridding
578 : void gwgrid(const Double*,
579 : Double*,
580 : const Complex*,
581 : Int*,
582 : Int*,
583 : Int*,
584 : const Int*,
585 : const Int*,
586 : const Float*,
587 : Int*,
588 : Int*,
589 : Double*,
590 : Double*,
591 : DComplex*,
592 : Int*,
593 : Int*,
594 : Int *,
595 : Int *,
596 : const Double*,
597 : const Double*,
598 : Int*,
599 : Int*,
600 : Int*,
601 : Int*,
602 : const Complex*,
603 : Int*,
604 : Int*,
605 : Double*);
606 :
607 : void sectgwgridd(const Double*,
608 : const Complex*,
609 : const Int*,
610 : const Int*,
611 : const Int*,
612 : const Int*,
613 : const Int*,
614 : const Float*,
615 : const Int*,
616 : DComplex*,
617 : const Int*,
618 : const Int*,
619 : const Int *,
620 : const Int *,
621 : //support
622 : const Int*,
623 : const Int*,
624 : const Int*,
625 : const Int*,
626 : const Complex*,
627 : const Int*,
628 : const Int*,
629 : Double*,
630 : //x0
631 : const Int*,
632 : const Int*,
633 : const Int*,
634 : const Int*,
635 : const Int*,
636 : const Int*,
637 : const Int*,
638 : const Int*,
639 : const Complex*
640 : );
641 :
642 : //Single precision gridding
643 : void sectgwgrids(const Double*,
644 : const Complex*,
645 : const Int*,
646 : const Int*,
647 : const Int*,
648 : const Int*,
649 : const Int*,
650 : const Float*,
651 : const Int*,
652 : Complex*,
653 : const Int*,
654 : const Int*,
655 : const Int *,
656 : const Int *,
657 : //support
658 : const Int*,
659 : const Int*,
660 : const Int*,
661 : const Int*,
662 : const Complex*,
663 : const Int*,
664 : const Int*,
665 : Double*,
666 : //x0
667 : const Int*,
668 : const Int*,
669 : const Int*,
670 : const Int*,
671 : const Int*,
672 : const Int*,
673 : const Int*,
674 : const Int*,
675 : const Complex*
676 : );
677 :
678 :
679 : void gwproj(const Double*,
680 : Double*,
681 : const Complex*,
682 : Int*,
683 : Int*,
684 : Int*,
685 : const Int*,
686 : const Int*,
687 : const Float*,
688 : Int*,
689 : Int*,
690 : Double*,
691 : Double*,
692 : Complex*,
693 : Int*,
694 : Int*,
695 : Int *,
696 : Int *,
697 : const Double*,
698 : const Double*,
699 : Int*,
700 : Int*,
701 : Int*,
702 : Int*,
703 : const Complex*,
704 : Int*,
705 : Int*,
706 : Double*);
707 :
708 : void sectdwgrid(const Double*,
709 : Complex*,
710 : const Int*,
711 : const Int*,
712 : const Int*,
713 : const Int*,
714 : const Int*,
715 : const Complex*,
716 : const Int*,
717 : const Int*,
718 : const Int *,
719 : const Int *,
720 : //support
721 : const Int*,
722 : const Int*,
723 : const Int*,
724 : const Int*,
725 : const Complex*,
726 : const Int*,
727 : const Int*,
728 : //rbeg, rend, loc, off, phasor
729 : const Int*,
730 : const Int*,
731 : const Int*,
732 : const Int*,
733 : const Complex*);
734 : void dwproj(const Double*,
735 : Double*,
736 : Complex*,
737 : Int*,
738 : Int*,
739 : const Int*,
740 : const Int*,
741 : Int*,
742 : Int*,
743 : Double*,
744 : Double*,
745 : const Complex*,
746 : Int*,
747 : Int*,
748 : Int *,
749 : Int *,
750 : const Double*,
751 : const Double*,
752 : Int*,
753 : Int*,
754 : Int*,
755 : Int*,
756 : const Complex*,
757 : Int*,
758 : Int*);
759 : }
760 5692 : void WProjectFT::put(const VisBuffer2& vb, Int row, Bool dopsf,
761 : FTMachine::Type type)
762 : {
763 :
764 :
765 : //Check if ms has changed then cache new spw and chan selection
766 : //if(vb.isNewMs())
767 : // matchAllSpwChans(vb);
768 :
769 : //Here we redo the match or use previous match
770 :
771 : //Channel matching for the actual spectral window of buffer
772 : //if(doConversion_p[vb.spectralWindows()[0]]){
773 5692 : matchChannel(vb);
774 : //}
775 : //else{
776 : // chanMap.resize();
777 : // chanMap=multiChanMap_p[vb.spectralWindows()[0]];
778 : //}
779 :
780 : //No point in reading data if its not matching in frequency
781 5692 : if(max(chanMap)==-1)
782 0 : return;
783 5692 : Timer tim;
784 5692 : tim.mark();
785 :
786 : //const Matrix<Float> *imagingweight;
787 : //imagingweight=&(vb.imagingWeight());
788 11384 : Matrix<Float> imagingweight;
789 5692 : getImagingWeight(imagingweight, vb);
790 5692 : if(dopsf) type=FTMachine::PSF;
791 :
792 11384 : Cube<Complex> data;
793 : //Fortran gridder need the flag as ints
794 11384 : Cube<Int> flags;
795 11384 : Matrix<Float> elWeight;
796 5692 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
797 :
798 :
799 : Bool iswgtCopy;
800 : const Float *wgtStorage;
801 5692 : wgtStorage=elWeight.getStorage(iswgtCopy);
802 :
803 :
804 : Bool isCopy;
805 5692 : const Complex *datStorage=0;
806 :
807 5692 : if(!dopsf)
808 3575 : datStorage=data.getStorage(isCopy);
809 :
810 :
811 : // If row is -1 then we pass through all rows
812 : Int startRow, endRow, nRow;
813 5692 : if (row==-1) {
814 5692 : nRow=vb.nRows();
815 5692 : startRow=0;
816 5692 : endRow=nRow-1;
817 : } else {
818 0 : nRow=1;
819 0 : startRow=row;
820 0 : endRow=row;
821 : }
822 :
823 : // Get the uvws in a form that Fortran can use and do that
824 : // necessary phase rotation. On a Pentium Pro 200 MHz
825 : // when null, this step takes about 50us per uvw point. This
826 : // is just barely noticeable for Stokes I continuum and
827 : // irrelevant for other cases.
828 11384 : Matrix<Double> uvw(negateUV(vb));
829 11384 : Vector<Double> dphase(vb.nRows());
830 5692 : dphase=0.0;
831 :
832 5692 : rotateUVW(uvw, dphase, vb);
833 5692 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
834 :
835 :
836 :
837 : // Take care of translation of Bools to Integer
838 5692 : Int idopsf=0;
839 5692 : if(dopsf) idopsf=1;
840 :
841 11384 : Vector<Int> rowFlags(vb.nRows());
842 5692 : rowFlags=0;
843 5692 : rowFlags(vb.flagRow())=true;
844 5692 : if(!usezero_p) {
845 2003584 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
846 1997892 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
847 : }
848 : }
849 :
850 :
851 : Bool del;
852 : // IPosition s(flags.shape());
853 11384 : Vector<Int> s(flags.shape().nelements());
854 5692 : convertArray(s, flags.shape().asVector());
855 5692 : Int nvp=s[0];
856 5692 : Int nvc=s[1];
857 5692 : Int nvisrow=s[2];
858 :
859 5692 : Int csamp=convSampling;
860 :
861 : Bool uvwcopy;
862 5692 : const Double *uvwstor=uvw.getStorage(uvwcopy);
863 : Bool gridcopy;
864 : Bool convcopy;
865 5692 : const Complex *convstor=convFunc.getStorage(convcopy);
866 11384 : Cube<Int> loc(3, nvc, nRow);
867 11384 : Cube<Int> off(3, nvc, nRow);
868 11384 : Matrix<Complex> phasor(nvc, nRow);
869 : Bool delphase;
870 5692 : Complex * phasorstor=phasor.getStorage(delphase);
871 5692 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
872 5692 : const Double * scalestor=uvScale.getStorage(del);
873 5692 : const Double * offsetstor=uvOffset.getStorage(del);
874 : Bool delloc, deloff;
875 5692 : Int * locstor=loc.getStorage(delloc);
876 5692 : Int * offstor=off.getStorage(deloff);
877 5692 : const Double *dpstor=dphase.getStorage(del);
878 5692 : Double cinv=Double(1.0)/C::c;
879 : Int irow;
880 5692 : Int dow=1;
881 5692 : Int nth=1;
882 : #ifdef _OPENMP
883 5692 : if(numthreads_p >0){
884 0 : nth=min(numthreads_p, omp_get_max_threads());
885 : }
886 : else{
887 5692 : nth= omp_get_max_threads();
888 : }
889 : //nth=min(4,nth);
890 : #endif
891 :
892 :
893 5692 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvc, scalestor, offsetstor, csamp, phasorstor, uvwstor, locstor, offstor, dpstor, cinv, dow) shared(startRow, endRow) num_threads(nth)
894 :
895 : {
896 : #pragma omp for schedule(dynamic)
897 : for (irow=startRow; irow<=endRow;irow++){
898 : //locateuvw(uvwstor,dpstor, visfreqstor, nvc, scalestor, offsetstor, csamp,
899 : // locstor,
900 : // offstor, phasorstor, irow, true);
901 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
902 : }
903 :
904 : }//end pragma parallel
905 5692 : Int maxx=0;
906 5692 : Int minx=1000000;
907 5692 : Int maxy=0;
908 5692 : Int miny=1000000;
909 : //Int maxsupp=max(convSupport);
910 5692 : loc.putStorage(locstor, delloc);
911 5692 : maxx=max(loc.yzPlane(0));
912 5692 : maxx=min(nx-1, maxx-1);
913 5692 : minx=min(loc.yzPlane(0));
914 5692 : minx=max(0, minx-1);
915 5692 : maxy=max(loc.yzPlane(1));
916 5692 : maxy=min(ny-1, maxy-1);
917 5692 : miny=min(loc.yzPlane(1));
918 5692 : miny=max(0,miny-1);
919 5692 : locstor=loc.getStorage(delloc);
920 : //cerr << "dels " << delloc << " " << deloff << endl;
921 : //cerr << "LOC " << min(loc.xzPlane(0)) << " max " << loc.shape() << endl;
922 5692 : timemass_p +=tim.real();
923 : Int ixsub, iysub, icounter;
924 5692 : ixsub=1;
925 5692 : iysub=1;
926 5692 : if (nth >4){
927 0 : ixsub=8;
928 0 : iysub=8;
929 : }
930 : else {
931 5692 : ixsub=2;
932 5692 : iysub=2;
933 : }
934 : //nxsub=nx;
935 : //nysub=ny;
936 5692 : Int rbeg=startRow+1;
937 5692 : Int rend=endRow+1;
938 :
939 5692 : const Int* pmapstor=polMap.getStorage(del);
940 5692 : const Int *cmapstor=chanMap.getStorage(del);
941 5692 : Int nc=nchan;
942 5692 : Int np=npol;
943 : // Int nxp=nx;
944 : // Int nyp=ny;
945 5692 : minx=0;
946 5692 : miny=0;
947 5692 : Int nxp=nx;
948 5692 : Int nyp=ny;
949 5692 : Int nxcopy=nx;
950 5692 : Int nycopy=ny;
951 :
952 5692 : Int csize=convSize;
953 5692 : Int wcsize=wConvSize;
954 5692 : const Int * flagstor=flags.getStorage(del);
955 5692 : const Int * rowflagstor=rowFlags.getStorage(del);
956 5692 : const Int * suppstor=convSupport.getStorage(del);
957 5692 : tim.mark();
958 11384 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
959 28460 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
960 22768 : sumwgt[icounter].resize(sumWeight.shape());
961 22768 : sumwgt[icounter].set(0.0);
962 : }
963 5692 : if(doneThreadPartition_p < 0){
964 47 : xsect_p.resize(ixsub*iysub);
965 47 : ysect_p.resize(ixsub*iysub);
966 47 : nxsect_p.resize(ixsub*iysub);
967 47 : nysect_p.resize(ixsub*iysub);
968 235 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
969 188 : findGridSector(nxp, nyp, ixsub, iysub, minx, miny, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
970 : }
971 : }
972 11384 : Vector<Int> xsect, ysect, nxsect, nysect;
973 5692 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
974 :
975 5692 : if(!useDoubleGrid_p){
976 0 : Complex *gridstor=griddedData.getStorage(gridcopy);
977 0 : #pragma omp parallel default(none) private(icounter, del) firstprivate(idopsf, uvwstor, datStorage, wgtStorage, flagstor, rowflagstor, convstor, pmapstor, cmapstor, gridstor, suppstor, nxp, nyp, nxcopy, nycopy, np, nc,ixsub, iysub, rend, rbeg, csamp, csize, wcsize, nvp, nvc, nvisrow, phasorstor, locstor, offstor, minx, miny, xsect, ysect, nxsect, nysect) shared(sumwgt) num_threads(nth)
978 : {
979 :
980 : #pragma omp for schedule(dynamic)
981 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
982 :
983 :
984 : Int x0=xsect(icounter);
985 : Int y0=ysect(icounter);
986 : Int nxsub=nxsect(icounter);
987 : Int nysub=nysect(icounter);
988 :
989 : sectgwgrids(uvwstor,
990 : datStorage,
991 : &nvp,
992 : &nvc,
993 : &idopsf,
994 : flagstor,
995 : rowflagstor,
996 : wgtStorage,
997 : &nvisrow,
998 : gridstor,
999 : &nxcopy,
1000 : &nycopy,
1001 : &np,
1002 : &nc,
1003 : suppstor,
1004 : &csize,
1005 : &csamp,
1006 : &wcsize,
1007 : convstor,
1008 : cmapstor,
1009 : pmapstor,
1010 : (sumwgt[icounter]).getStorage(del),
1011 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
1012 : phasorstor);
1013 : }
1014 : }//end pragma parallel
1015 0 : if(dopsf && (nth > 4))
1016 0 : tweakGridSector(nx, ny, ixsub, iysub);
1017 0 : timegrid_p+=tim.real();
1018 :
1019 0 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1020 0 : sumWeight=sumWeight+sumwgt[icounter];
1021 : }
1022 0 : griddedData.putStorage(gridstor, gridcopy);
1023 :
1024 : }
1025 : else{
1026 5692 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
1027 5692 : #pragma omp parallel default(none) private(icounter,del) firstprivate(idopsf, uvwstor, datStorage, wgtStorage, flagstor, rowflagstor, convstor, pmapstor, cmapstor, gridstor, suppstor, nxp, nyp, nxcopy, nycopy, np, nc,ixsub, iysub, rend, rbeg, csamp, csize, wcsize, nvp, nvc, nvisrow, phasorstor, locstor, offstor,minx,miny, xsect, ysect, nxsect, nysect) shared(sumwgt) num_threads(nth)
1028 : {
1029 : #pragma omp for schedule(dynamic)
1030 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
1031 :
1032 : Int x0=xsect(icounter);
1033 : Int y0=ysect(icounter);
1034 : Int nxsub=nxsect(icounter);
1035 : Int nysub=nysect(icounter);
1036 :
1037 : sectgwgridd(uvwstor,
1038 : datStorage,
1039 : &nvp,
1040 : &nvc,
1041 : &idopsf,
1042 : flagstor,
1043 : rowflagstor,
1044 : wgtStorage,
1045 : &nvisrow,
1046 : gridstor,
1047 : &nxcopy,
1048 : &nycopy,
1049 : &np,
1050 : &nc,
1051 : suppstor,
1052 : &csize,
1053 : &csamp,
1054 : &wcsize,
1055 : convstor,
1056 : cmapstor,
1057 : pmapstor,
1058 : (sumwgt[icounter]).getStorage(del),
1059 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
1060 : phasorstor);
1061 : }
1062 : }//end pragma parallel
1063 5692 : if(dopsf && (nth > 4))
1064 0 : tweakGridSector(nx, ny, ixsub, iysub);
1065 5692 : timegrid_p+=tim.real();
1066 :
1067 28460 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1068 22768 : sumWeight=sumWeight+sumwgt[icounter];
1069 : }
1070 5692 : griddedData2.putStorage(gridstor, gridcopy);
1071 : }
1072 5692 : uvw.freeStorage(uvwstor, uvwcopy);
1073 5692 : convFunc.freeStorage(convstor, convcopy);
1074 :
1075 5692 : if(!dopsf)
1076 3575 : data.freeStorage(datStorage, isCopy);
1077 5692 : elWeight.freeStorage(wgtStorage,iswgtCopy);
1078 : }
1079 :
1080 :
1081 :
1082 1458 : void WProjectFT::get(VisBuffer2& vb, Int row)
1083 : {
1084 :
1085 :
1086 1458 : findConvFunction(*image, vb);
1087 : // If row is -1 then we pass through all rows
1088 : Int startRow, endRow, nRow;
1089 1458 : if (row==-1) {
1090 1458 : nRow=vb.nRows();
1091 1458 : startRow=0;
1092 1458 : endRow=nRow-1;
1093 : //vb.modelVisCube()=Complex(0.0,0.0);
1094 : } else {
1095 0 : nRow=1;
1096 0 : startRow=row;
1097 0 : endRow=row;
1098 : //vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
1099 : }
1100 :
1101 :
1102 : //Channel matching for the actual spectral window of buffer
1103 1458 : matchChannel(vb);
1104 :
1105 : //No point in reading data if its not matching in frequency
1106 1458 : if(max(chanMap)==-1)
1107 0 : return;
1108 :
1109 :
1110 : // Get the uvws in a form that Fortran can use
1111 2916 : Matrix<Double> uvw(negateUV(vb));
1112 2916 : Vector<Double> dphase(vb.nRows());
1113 1458 : dphase=0.0;
1114 :
1115 1458 : rotateUVW(uvw, dphase, vb);
1116 1458 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1117 :
1118 : // This is the convention for dphase
1119 : // dphase*=-1.0;
1120 :
1121 :
1122 : //Check if ms has changed then cache new spw and chan selection
1123 : //if(vb.isNewMs())
1124 : // matchAllSpwChans(vb);
1125 : //Here we redo the match or use previous match
1126 :
1127 :
1128 2916 : Cube<Complex> data;
1129 2916 : Cube<Int> flags;
1130 1458 : getInterpolateArrays(vb, data, flags);
1131 :
1132 :
1133 :
1134 : Complex *datStorage;
1135 : Bool isCopy;
1136 1458 : datStorage=data.getStorage(isCopy);
1137 :
1138 2916 : Vector<Int> rowFlags(vb.nRows());
1139 1458 : rowFlags=0;
1140 1458 : rowFlags(vb.flagRow())=true;
1141 1458 : if(!usezero_p) {
1142 513216 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1143 511758 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
1144 : }
1145 : }
1146 :
1147 1458 : Int nvp=data.shape()(0);
1148 1458 : Int nvc=data.shape()(1);
1149 1458 : Int nvisrow=data.shape()(2);
1150 1458 : Int nc=nchan;
1151 1458 : Int np=npol;
1152 1458 : Int nxp=nx;
1153 1458 : Int nyp=ny;
1154 2916 : Cube<Int> loc(3, nvc, nvisrow);
1155 2916 : Cube<Int> off(3, nvc, nRow);
1156 1458 : Int csamp=convSampling;
1157 1458 : Int csize=convSize;
1158 1458 : Int wcsize=wConvSize;
1159 2916 : Matrix<Complex> phasor(nvc, nRow);
1160 : Bool delphase;
1161 : Bool del;
1162 1458 : Complex * phasorstor=phasor.getStorage(delphase);
1163 1458 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1164 1458 : const Double * scalestor=uvScale.getStorage(del);
1165 1458 : const Double * offsetstor=uvOffset.getStorage(del);
1166 1458 : Int * locstor=loc.getStorage(del);
1167 1458 : Int * offstor=off.getStorage(del);
1168 1458 : const Int * flagstor=flags.getStorage(del);
1169 1458 : const Int * rowflagstor=rowFlags.getStorage(del);
1170 1458 : const Double *dpstor=dphase.getStorage(del);
1171 : Bool uvwcopy;
1172 1458 : const Double *uvwstor=uvw.getStorage(uvwcopy);
1173 : Bool gridcopy;
1174 1458 : const Complex *gridstor=griddedData.getStorage(gridcopy);
1175 : Bool convcopy;
1176 1458 : const Complex *convstor=convFunc.getStorage(convcopy);
1177 1458 : const Int* pmapstor=polMap.getStorage(del);
1178 1458 : const Int *cmapstor=chanMap.getStorage(del);
1179 1458 : const Int * suppstor=convSupport.getStorage(del);
1180 : Int irow;
1181 1458 : Int nth=1;
1182 : #ifdef _OPENMP
1183 1458 : if(numthreads_p >0){
1184 0 : nth=min(numthreads_p, omp_get_max_threads());
1185 : }
1186 : else{
1187 1458 : nth= omp_get_max_threads();
1188 : }
1189 : // nth=min(4,nth);
1190 : #endif
1191 1458 : Int dow=1;
1192 1458 : Double cinv=Double(1.0)/C::c;
1193 1458 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvc, scalestor, offsetstor, csamp, phasorstor, uvwstor, locstor, offstor, dpstor, dow, cinv) shared(startRow, endRow) num_threads(nth)
1194 : {
1195 :
1196 : #pragma omp for schedule(dynamic)
1197 : for (irow=startRow; irow<=endRow; ++irow){
1198 : /*locateuvw(uvwstor,dpstor, visfreqstor, nvc, scalestor, offsetstor, csamp,
1199 : locstor,
1200 : offstor, phasorstor, irow, true);*/
1201 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
1202 : }
1203 :
1204 : }//end pragma parallel
1205 1458 : Int rbeg=startRow+1;
1206 1458 : Int rend=endRow+1;
1207 1458 : Int npart=nth;
1208 1458 : Timer tim;
1209 1458 : tim.mark();
1210 :
1211 1458 : Int ix=0;
1212 1458 : #pragma omp parallel default(none) private(ix, rbeg, rend) firstprivate(uvwstor, datStorage, flagstor, rowflagstor, convstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc, suppstor, csamp, csize, wcsize, nvp, nvc, nvisrow, phasorstor, locstor, offstor) shared(npart) num_threads(nth)
1213 : {
1214 :
1215 : #pragma omp for schedule(dynamic)
1216 : for (ix=0; ix< npart; ++ix){
1217 : rbeg=ix*(nvisrow/npart)+1;
1218 : rend=(ix != (npart-1)) ? (rbeg+(nvisrow/npart)-1) : (rbeg+(nvisrow/npart)+nvisrow%npart-1) ;
1219 : // cerr << "rbeg " << rbeg << " rend " << rend << " nRow " << nvisrow << endl;
1220 : sectdwgrid(uvwstor,
1221 : datStorage,
1222 : &nvp,
1223 : &nvc,
1224 : flagstor,
1225 : rowflagstor,
1226 : &nvisrow,
1227 : gridstor,
1228 : &nxp,
1229 : &nyp,
1230 : &np,
1231 : &nc,
1232 : suppstor,
1233 : &csize,
1234 : &csamp,
1235 : &wcsize,
1236 : convstor,
1237 : cmapstor,
1238 : pmapstor,
1239 : &rbeg, &rend, locstor, offstor, phasorstor);
1240 : }
1241 :
1242 : }//end pragma parallel
1243 1458 : data.putStorage(datStorage, isCopy);
1244 1458 : uvw.freeStorage(uvwstor, uvwcopy);
1245 1458 : griddedData.freeStorage(gridstor, gridcopy);
1246 1458 : convFunc.freeStorage(convstor, convcopy);
1247 1458 : timedegrid_p+=tim.real();
1248 :
1249 1458 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1250 : }
1251 :
1252 :
1253 :
1254 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1255 : // return the resulting image
1256 73 : ImageInterface<Complex>& WProjectFT::getImage(Matrix<Float>& weights,
1257 : Bool normalize)
1258 : {
1259 : //AlwaysAssert(lattice, AipsError);
1260 73 : AlwaysAssert(image, AipsError);
1261 :
1262 73 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
1263 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
1264 : }
1265 73 : logIO() << LogOrigin("WProjectFT", "getImage") << LogIO::NORMAL;
1266 :
1267 73 : weights.resize(sumWeight.shape());
1268 :
1269 73 : convertArray(weights, sumWeight);
1270 :
1271 : // If the weights are all zero then we cannot normalize
1272 : // otherwise we don't care.
1273 73 : if(max(weights)==0.0) {
1274 0 : if(normalize) {
1275 : logIO() << LogIO::SEVERE
1276 : << "No useful data in WProjectFT: weights all zero"
1277 0 : << LogIO::POST;
1278 : }
1279 : else {
1280 : logIO() << LogIO::WARN
1281 : << "No useful data in WProjectFT: weights all zero"
1282 0 : << LogIO::POST;
1283 : }
1284 : }
1285 : else {
1286 : logIO() << LogIO::DEBUGGING
1287 73 : << "Starting FFT and scaling of image" << LogIO::POST;
1288 :
1289 73 : if(useDoubleGrid_p){
1290 146 : ArrayLattice<DComplex> darrayLattice(griddedData2);
1291 73 : LatticeFFT::cfft2d(darrayLattice,false);
1292 73 : griddedData.resize(griddedData2.shape());
1293 73 : convertArray(griddedData, griddedData2);
1294 73 : SynthesisUtilMethods::getResource("mem peak in getImage");
1295 73 : griddedData2.resize();
1296 73 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1297 73 : lattice=arrayLattice;
1298 : }else{
1299 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1300 0 : lattice=arrayLattice;
1301 0 : LatticeFFT::cfft2d(*lattice,false);
1302 :
1303 : }
1304 :
1305 :
1306 146 : const IPosition latticeShape = lattice->shape();
1307 :
1308 :
1309 : {
1310 73 : Int inx = lattice->shape()(0);
1311 73 : Int iny = lattice->shape()(1);
1312 146 : Vector<Complex> correction(inx);
1313 73 : correction=Complex(1.0, 0.0);
1314 :
1315 73 : Int npixCorr= max(nx,ny);
1316 146 : Vector<Float> sincConv(npixCorr);
1317 78113 : for (Int ix=0;ix<npixCorr;ix++) {
1318 78040 : Float x=C::pi*Float(ix-npixCorr/2)/(Float(npixCorr)*Float(convSampling));
1319 78040 : if(ix==npixCorr/2) {
1320 73 : sincConv(ix)=1.0;
1321 : }
1322 : else {
1323 77967 : sincConv(ix)=sin(x)/x;
1324 : }
1325 : }
1326 :
1327 146 : IPosition cursorShape(4, inx, 1, 1, 1);
1328 146 : IPosition axisPath(4, 0, 1, 2, 3);
1329 146 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1330 146 : LatticeIterator<Complex> lix(*lattice, lsx);
1331 109513 : for(lix.reset();!lix.atEnd();lix++) {
1332 109440 : Int pol=lix.position()(2);
1333 109440 : Int chan=lix.position()(3);
1334 109440 : if(weights(pol, chan)!=0.0) {
1335 109440 : Int iy=lix.position()(1);
1336 109440 : gridder->correctX1D(correction,iy);
1337 155628840 : for (Int ix=0;ix<nx;ix++) {
1338 155519400 : correction(ix)*=sincConv(ix)*sincConv(iy);
1339 : }
1340 109440 : lix.rwVectorCursor()/=correction;
1341 109440 : if(normalize) {
1342 0 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
1343 0 : lix.rwCursor()*=rnorm;
1344 : }
1345 : else {
1346 109440 : Complex rnorm(Float(inx)*Float(iny));
1347 109440 : lix.rwCursor()*=rnorm;
1348 : }
1349 : }
1350 : else {
1351 0 : lix.woCursor()=0.0;
1352 : }
1353 : }
1354 : }
1355 :
1356 73 : if(!isTiled) {
1357 146 : LatticeLocker lock1 (*(image), FileLocker::Write);
1358 : // Check the section from the image BEFORE converting to a lattice
1359 146 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
1360 219 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1361 146 : IPosition stride(4, 1);
1362 146 : IPosition trc(blc+image->shape()-stride);
1363 :
1364 : // Do the copy
1365 73 : image->put(griddedData(blc, trc));
1366 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
1367 73 : griddedData.resize(IPosition(1,0));
1368 : }
1369 : }
1370 73 : if(!lattice.null()) lattice=nullptr;
1371 73 : if(!arrayLattice.null()) lattice=nullptr;
1372 73 : griddedData.resize();
1373 73 : return *image;
1374 : }
1375 :
1376 : // Get weight image
1377 0 : void WProjectFT::getWeightImage(ImageInterface<Float>& weightImage,
1378 : Matrix<Float>& weights)
1379 : {
1380 :
1381 0 : logIO() << LogOrigin("WProjectFT", "getWeightImage") << LogIO::NORMAL;
1382 :
1383 0 : weights.resize(sumWeight.shape());
1384 0 : convertArray(weights,sumWeight);
1385 :
1386 0 : const IPosition latticeShape = weightImage.shape();
1387 :
1388 0 : Int inx=latticeShape(0);
1389 0 : Int iny=latticeShape(1);
1390 :
1391 0 : IPosition loc(2, 0);
1392 0 : IPosition cursorShape(4, inx, iny, 1, 1);
1393 0 : IPosition axisPath(4, 0, 1, 2, 3);
1394 0 : LatticeStepper lsx(latticeShape, cursorShape, axisPath);
1395 0 : LatticeIterator<Float> lix(weightImage, lsx);
1396 0 : for(lix.reset();!lix.atEnd();lix++) {
1397 0 : Int pol=lix.position()(2);
1398 0 : Int chan=lix.position()(3);
1399 0 : lix.rwCursor()=weights(pol,chan);
1400 : }
1401 0 : }
1402 :
1403 0 : Bool WProjectFT::toRecord(String& error,
1404 : RecordInterface& outRec, Bool withImage, const String diskimage)
1405 : {
1406 :
1407 : /*
1408 :
1409 : CountedPtr<WPConvFunc> wpConvFunc_p;
1410 : */
1411 :
1412 : // Save the current WProjectFT object to an output state record
1413 0 : Bool retval = true;
1414 : //save the base class variables
1415 : //this is a memory hog and slow on saving and recovering...better to recompute convfunctions
1416 : /* Record wpconvrec;
1417 : if(wpConvFunc_p->toRecord(wpconvrec))
1418 : outRec.defineRecord("wpconvfunc", wpconvrec);
1419 : */
1420 0 : Float elpadd=padding_p;
1421 0 : if(toVis_p && withImage)
1422 0 : elpadd=1.0;
1423 0 : if(!FTMachine::toRecord(error, outRec, withImage, diskimage))
1424 0 : return false;
1425 :
1426 0 : outRec.define("uvscale", uvScale);
1427 0 : outRec.define("uvoffset", uvOffset);
1428 0 : outRec.define("istiled", isTiled);
1429 0 : outRec.define("cachesize", Int64(cachesize));
1430 0 : outRec.define("tilesize", tilesize);
1431 0 : outRec.define("maxabsdata", maxAbsData);
1432 0 : Vector<Int> center_loc(4), offset_loc(4);
1433 0 : for (Int k=0; k<4 ; k++){
1434 0 : center_loc(k)=centerLoc(k);
1435 0 : offset_loc(k)=offsetLoc(k);
1436 : }
1437 0 : outRec.define("centerloc", center_loc);
1438 0 : outRec.define("offsetloc", offset_loc);
1439 0 : outRec.define("padding", elpadd);
1440 0 : outRec.define("nwplanes", nWPlanes_p);
1441 0 : outRec.define("savedwscale", savedWScale_p);
1442 0 : outRec.define("usezero", usezero_p);
1443 : ///no need to save convfunc as it can be big and is recalculated anyways
1444 : ///outRec.define("convfunc", convFunc);
1445 0 : outRec.define("convsampling", convSampling);
1446 0 : outRec.define("convsize", convSize);
1447 0 : outRec.define("convsupport", convSupport);
1448 0 : outRec.define("convsizes", convSizes_p);
1449 0 : outRec.define("wconvsize", wConvSize);
1450 0 : outRec.define("lastindex", lastIndex_p);
1451 0 : outRec.define("minw", minW_p);
1452 0 : outRec.define("maxw", maxW_p);
1453 0 : outRec.define("rmsw", rmsW_p);
1454 :
1455 :
1456 :
1457 0 : return retval;
1458 : }
1459 :
1460 0 : Bool WProjectFT::fromRecord(String& error,
1461 : const RecordInterface& inRec)
1462 : {
1463 0 : if(!FTMachine::fromRecord(error, inRec))
1464 0 : return false;
1465 0 : machineName_p="WProjectFT";
1466 0 : Bool retval = true;
1467 0 : imageCache=0; lattice=0; arrayLattice=0;
1468 0 : inRec.get("uvscale", uvScale);
1469 0 : inRec.get("uvoffset", uvOffset);
1470 0 : inRec.get("istiled", isTiled);
1471 0 : cachesize=inRec.asInt64("cachesize");
1472 0 : inRec.get("tilesize", tilesize);
1473 0 : inRec.get("maxabsdata", maxAbsData);
1474 0 : Vector<Int> center_loc(4), offset_loc(4);
1475 0 : inRec.get("centerloc", center_loc);
1476 0 : inRec.get("offsetloc", offset_loc);
1477 0 : uInt ndim4 = 4;
1478 0 : centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2),
1479 0 : center_loc(3));
1480 0 : offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2),
1481 0 : offset_loc(3));
1482 0 : inRec.get("minw", minW_p);
1483 0 : inRec.get("maxw", maxW_p);
1484 0 : inRec.get("rmsw", rmsW_p);
1485 0 : if(inRec.isDefined("wpconvfunc")){
1486 0 : wpConvFunc_p=new WPConvFunc(inRec.asRecord("wpconvfunc"));
1487 : }
1488 : else{
1489 0 : wpConvFunc_p=new WPConvFunc(minW_p, maxW_p, rmsW_p);
1490 : }
1491 0 : inRec.get("padding", padding_p);
1492 0 : inRec.get("nwplanes", nWPlanes_p);
1493 0 : inRec.get("savedwscale", savedWScale_p);
1494 0 : inRec.get("usezero", usezero_p);
1495 : //inRec.get("convfunc", convFunc);
1496 0 : convFunc.resize();
1497 0 : inRec.get("convsampling", convSampling);
1498 0 : inRec.get("convsize", convSize);
1499 0 : inRec.get("convsupport", convSupport);
1500 0 : inRec.get("convsizes", convSizes_p);
1501 0 : inRec.get("wconvsize", wConvSize);
1502 0 : inRec.get("lastindex", lastIndex_p);
1503 0 : gridder=0;
1504 : ///setup some of the parameters
1505 0 : init();
1506 :
1507 :
1508 :
1509 0 : return retval;
1510 : }
1511 :
1512 14313 : void WProjectFT::ok() {
1513 14313 : AlwaysAssert(image, AipsError);
1514 14313 : }
1515 :
1516 : // Make a plain straightforward honest-to-God image. This returns
1517 : // a complex image, without conversion to Stokes. The representation
1518 : // is that required for the visibilities.
1519 : //----------------------------------------------------------------------
1520 0 : void WProjectFT::makeImage(FTMachine::Type type,
1521 : VisibilityIterator2& vi,
1522 : ImageInterface<Complex>& theImage,
1523 : Matrix<Float>& weight) {
1524 :
1525 :
1526 0 : logIO() << LogOrigin("WProjectFT", "makeImage") << LogIO::NORMAL;
1527 :
1528 0 : if(type==FTMachine::COVERAGE) {
1529 : logIO() << "Type COVERAGE not defined for Fourier transforms"
1530 0 : << LogIO::EXCEPTION;
1531 : }
1532 :
1533 :
1534 :
1535 : // Loop over all visibilities and pixels
1536 0 : VisBuffer2 *vb=vi.getVisBuffer();
1537 :
1538 : // Initialize put (i.e. transform to Sky) for this model
1539 0 : vi.origin();
1540 :
1541 0 : if(vb->polarizationFrame()==MSIter::Linear) {
1542 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
1543 : }
1544 : else {
1545 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
1546 : }
1547 :
1548 0 : initializeToSky(theImage,weight,*vb);
1549 :
1550 : // Loop over the visibilities, putting VisBuffers
1551 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1552 0 : for (vi.origin(); vi.more(); vi.next()) {
1553 :
1554 0 : switch(type) {
1555 0 : case FTMachine::RESIDUAL:
1556 : //vb.visCube()=vb.correctedVisCube();
1557 0 : vb->setVisCube(vb->visCubeCorrected()-vb->visCubeModel());
1558 0 : put(*vb, -1, false);
1559 0 : break;
1560 0 : case FTMachine::MODEL:
1561 0 : vb->setVisCube(vb->visCubeModel());
1562 0 : put(*vb, -1, false);
1563 0 : break;
1564 0 : case FTMachine::CORRECTED:
1565 0 : vb->setVisCube(vb->visCubeCorrected());
1566 0 : put(*vb, -1, false);
1567 0 : break;
1568 0 : case FTMachine::PSF:
1569 0 : vb->setVisCube(Complex(1.0,0.0));
1570 0 : put(*vb, -1, true);
1571 0 : break;
1572 0 : case FTMachine::OBSERVED:
1573 : default:
1574 0 : put(*vb, -1, false);
1575 0 : break;
1576 : }
1577 : }
1578 : }
1579 0 : finalizeToSky();
1580 : // Normalize by dividing out weights, etc.
1581 0 : getImage(weight, true);
1582 0 : }
1583 :
1584 :
1585 127 : String WProjectFT::name() const {
1586 :
1587 127 : return machineName_p;
1588 :
1589 : }
1590 7 : void WProjectFT::wStat(vi::VisibilityIterator2& vi, Double& minW, Double& maxW, Double& rmsW){
1591 7 : VisBuffer2* vb= vi.getVisBuffer();
1592 7 : maxW=0.0;
1593 7 : minW=1e99;
1594 7 : Double nval=0;
1595 7 : rmsW=0.0;
1596 28 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
1597 273 : for (vi.origin();vi.more();vi.next()) {
1598 252 : maxW=max(maxW, max(abs(vb->uvw().row(2)*max(vb->getFrequencies(0))))/C::c);
1599 252 : minW=min(minW, min(abs(vb->uvw().row(2)*min(vb->getFrequencies(0))))/C::c);
1600 : ///////////some shenanigans as some compilers is confusing * operator for vector
1601 252 : Vector<Double> elw;
1602 252 : elw=vb->uvw().row(2);
1603 252 : elw*=vb->uvw().row(2);
1604 : //////////////////////////////////////////////////
1605 252 : rmsW+=sum(elw);
1606 :
1607 252 : nval+=Double(vb->nRows());
1608 : }
1609 : }
1610 7 : rmsW=sqrt(rmsW/Double(nval))*max(vb->getFrequencies(0))/C::c;
1611 :
1612 7 : }
1613 :
1614 :
1615 : } // end of namespace refim
1616 : } //# NAMESPACE CASA - END
1617 :
|
