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: casa-feedback@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 0 : : FTMachine(), padding_p(1.0), nWPlanes_p(nWPlanes),
94 0 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
95 0 : gridder(0), isTiled(false),
96 0 : 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 89 : WProjectFT::WProjectFT(Int nWPlanes,
107 : MPosition mLocation,
108 : Long icachesize, Int itilesize,
109 89 : Bool usezero, Float padding, Bool useDoublePrec, const Double minW, const Double maxW, const Double rmsW)
110 89 : : FTMachine(), padding_p(padding), nWPlanes_p(nWPlanes),
111 89 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
112 89 : gridder(0), isTiled(false),
113 89 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
114 89 : usezero_p(usezero),
115 178 : 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 89 : convSize=0;
119 89 : savedWScale_p=0.0;
120 89 : tangentSpecified_p=false;
121 89 : mLocation_p=mLocation;
122 89 : lastIndex_p=0;
123 89 : wpConvFunc_p=new WPConvFunc(minW, maxW, rmsW);
124 89 : useDoubleGrid_p=useDoublePrec;
125 89 : }
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 2 : : FTMachine(), padding_p(padding), nWPlanes_p(nWPlanes),
132 2 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
133 2 : gridder(0), isTiled(false),
134 2 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
135 2 : usezero_p(usezero),
136 4 : 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 94 : 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 94 : isTiled=false;
226 94 : CompositeNumber cn(uInt(image->shape()(0)*2));
227 94 : nx = cn.nextLargerEven(Int(padding_p*Float(image->shape()(0))-0.5));
228 94 : ny = cn.nextLargerEven(Int(padding_p*Float(image->shape()(1))-0.5));
229 94 : npol = image->shape()(2);
230 94 : 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 94 : isTiled=false;
241 :
242 :
243 94 : sumWeight.resize(npol, nchan);
244 :
245 94 : wConvSize=max(0, nWPlanes_p);
246 94 : convSupport.resize(wConvSize);
247 94 : convSupport=0;
248 :
249 94 : uvScale.resize(3);
250 94 : uvScale=0.0;
251 94 : uvScale(0)=Float(nx)*image->coordinates().increment()(0);
252 94 : uvScale(1)=Float(ny)*image->coordinates().increment()(1);
253 94 : if(savedWScale_p==0.0){
254 54 : uvScale(2)=Float(wConvSize)*abs(image->coordinates().increment()(0));
255 : }
256 : else{
257 40 : uvScale(2)=savedWScale_p;
258 : }
259 94 : uvOffset.resize(3);
260 94 : uvOffset(0)=nx/2;
261 94 : uvOffset(1)=ny/2;
262 94 : uvOffset(2)=0;
263 :
264 :
265 :
266 94 : if(gridder) delete gridder; gridder=0;
267 188 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
268 94 : uvScale, uvOffset,
269 94 : "SF");
270 :
271 : // Set up image cache needed for gridding.
272 94 : 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 94 : 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 0 : }
290 94 : }
291 :
292 : // This is nasty, we should use CountedPointers here.
293 240 : WProjectFT::~WProjectFT() {
294 120 : if(imageCache) delete imageCache; imageCache=0;
295 120 : 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 240 : }
310 :
311 :
312 45 : void WProjectFT::setConvFunc(CountedPtr<WPConvFunc>& pbconvFunc){
313 :
314 45 : wpConvFunc_p=pbconvFunc;
315 45 : }
316 45 : CountedPtr<WPConvFunc>& WProjectFT::getConvFunc(){
317 :
318 45 : return wpConvFunc_p;
319 : }
320 :
321 1824 : void WProjectFT::findConvFunction(const ImageInterface<Complex>& image,
322 : const VisBuffer2& vb) {
323 :
324 :
325 :
326 :
327 :
328 1824 : wpConvFunc_p->findConvFunction(image, vb, wConvSize, uvScale, uvOffset,
329 1824 : padding_p,
330 1824 : convSampling,
331 1824 : convFunc, convSize, convSupport,
332 1824 : savedWScale_p);
333 :
334 1824 : nWPlanes_p=convSupport.nelements();
335 1824 : wConvSize=max(1,nWPlanes_p);
336 1824 : uvScale(2)=savedWScale_p;
337 :
338 :
339 :
340 1824 : }
341 :
342 16 : void WProjectFT::initializeToVis(ImageInterface<Complex>& iimage,
343 : const VisBuffer2& vb)
344 : {
345 16 : image=&iimage;
346 16 : toVis_p=true;
347 16 : ok();
348 :
349 : // if(convSize==0) {
350 16 : init();
351 : // }
352 16 : findConvFunction(*image, vb);
353 :
354 :
355 : // Initialize the maps for polarization and channel. These maps
356 : // translate visibility indices into image indices
357 16 : 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 16 : 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 16 : prepGridForDegrid();
377 16 : }
378 :
379 16 : void WProjectFT::prepGridForDegrid(){
380 16 : IPosition gridShape(4, nx, ny, npol, nchan);
381 16 : griddedData.resize(gridShape);
382 16 : griddedData=Complex(0.0);
383 :
384 :
385 16 : IPosition stride(4, 1);
386 32 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
387 32 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
388 32 : IPosition trc(blc+image->shape()-stride);
389 :
390 16 : IPosition start(4, 0);
391 16 : griddedData(blc, trc) = image->getSlice(start, image->shape());
392 :
393 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
394 16 : arrayLattice = new ArrayLattice<Complex>(griddedData);
395 16 : lattice=arrayLattice;
396 :
397 16 : logIO() << LogIO::DEBUGGING << "Starting FFT of image" << LogIO::POST;
398 :
399 :
400 16 : Vector<Float> sincConvX(nx);
401 17556 : for (Int ix=0;ix<nx;ix++) {
402 17540 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
403 17540 : if(ix==nx/2) {
404 16 : sincConvX(ix)=1.0;
405 : }
406 : else {
407 17524 : sincConvX(ix)=sin(x)/x;
408 : }
409 : }
410 16 : Vector<Float> sincConvY(ny);
411 17556 : for (Int ix=0;ix<ny;ix++) {
412 17540 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
413 17540 : if(ix==ny/2) {
414 16 : sincConvY(ix)=1.0;
415 : }
416 : else {
417 17524 : sincConvY(ix)=sin(x)/x;
418 : }
419 : }
420 :
421 16 : Vector<Complex> correction(nx);
422 16 : correction=Complex(1.0, 0.0);
423 : // Do the Grid-correction
424 16 : IPosition cursorShape(4, nx, 1, 1, 1);
425 16 : IPosition axisPath(4, 0, 1, 2, 3);
426 16 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
427 16 : LatticeIterator<Complex> lix(*lattice, lsx);
428 26396 : for(lix.reset();!lix.atEnd();lix++) {
429 26380 : Int iy=lix.position()(1);
430 26380 : gridder->correctX1D(correction,iy);
431 36722180 : for (Int ix=0;ix<nx;ix++) {
432 36695800 : correction(ix)/=(sincConvX(ix)*sincConvY(iy));
433 : }
434 26380 : lix.rwVectorCursor()/=correction;
435 : }
436 :
437 : // Now do the FFT2D in place
438 16 : LatticeFFT::cfft2d(*lattice);
439 :
440 16 : logIO() << LogIO::DEBUGGING << "Finished FFT" << LogIO::POST;
441 :
442 :
443 :
444 16 : }
445 :
446 :
447 15 : void WProjectFT::finalizeToVis()
448 : {
449 15 : logIO() << LogOrigin("WProjectFT", "finalizeToVis") << LogIO::NORMAL;
450 15 : logIO() <<LogIO::NORMAL2<< "Time to degrid " << timedegrid_p << LogIO::POST ;
451 15 : timedegrid_p=0.0;
452 15 : if(!arrayLattice.null()) arrayLattice=0;
453 15 : if(!lattice.null()) lattice=0;
454 15 : griddedData.resize();
455 15 : 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 0 : }
466 15 : }
467 :
468 :
469 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
470 : // grid.
471 78 : void WProjectFT::initializeToSky(ImageInterface<Complex>& iimage,
472 : Matrix<Float>& weight,
473 : const VisBuffer2& vb)
474 : {
475 : // image always points to the image
476 78 : image=&iimage;
477 78 : toVis_p=false;
478 :
479 : // if(convSize==0) {
480 78 : init();
481 : // }
482 78 : findConvFunction(*image, vb);
483 :
484 :
485 : // Initialize the maps for polarization and channel. These maps
486 : // translate visibility indices into image indices
487 78 : 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 78 : isTiled=false;
501 78 : sumWeight=0.0;
502 78 : weight.resize(sumWeight.shape());
503 78 : 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 78 : 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 78 : IPosition gridShape(4, nx, ny, npol, nchan);
515 78 : if(!useDoubleGrid_p){
516 1 : griddedData.resize(gridShape);
517 1 : griddedData=Complex(0.0);
518 : }
519 : else{
520 : //griddedData.resize();
521 77 : griddedData2.resize(gridShape);
522 77 : griddedData2=DComplex(0.0);
523 : }
524 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
525 :
526 78 : }
527 : //AlwaysAssert(lattice, AipsError);
528 :
529 78 : }
530 :
531 78 : void WProjectFT::finalizeToSky()
532 : {
533 78 : logIO() << LogOrigin("WProjectFT", "finalizeToSky") << LogIO::NORMAL;
534 78 : logIO() <<LogIO::NORMAL2<< "Time to massage data " << timemass_p << LogIO::POST;
535 78 : logIO() << LogIO::NORMAL2 <<"Time to grid data " << timegrid_p << LogIO::POST;
536 78 : timemass_p=0.0;
537 78 : 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 78 : 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 0 : }
550 78 : }
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 6036 : 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 6036 : 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 6036 : if(max(chanMap)==-1)
782 0 : return;
783 6036 : Timer tim;
784 6036 : tim.mark();
785 :
786 : //const Matrix<Float> *imagingweight;
787 : //imagingweight=&(vb.imagingWeight());
788 6036 : Matrix<Float> imagingweight;
789 6036 : getImagingWeight(imagingweight, vb);
790 6036 : if(dopsf) type=FTMachine::PSF;
791 :
792 6036 : Cube<Complex> data;
793 : //Fortran gridder need the flag as ints
794 6036 : Cube<Int> flags;
795 6036 : Matrix<Float> elWeight;
796 6036 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
797 :
798 :
799 : Bool iswgtCopy;
800 : const Float *wgtStorage;
801 6036 : wgtStorage=elWeight.getStorage(iswgtCopy);
802 :
803 :
804 : Bool isCopy;
805 6036 : const Complex *datStorage=0;
806 :
807 6036 : if(!dopsf)
808 3883 : datStorage=data.getStorage(isCopy);
809 :
810 :
811 : // If row is -1 then we pass through all rows
812 : Int startRow, endRow, nRow;
813 6036 : if (row==-1) {
814 6036 : nRow=vb.nRows();
815 6036 : startRow=0;
816 6036 : 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 6036 : Matrix<Double> uvw(negateUV(vb));
829 6036 : Vector<Double> dphase(vb.nRows());
830 6036 : dphase=0.0;
831 :
832 6036 : rotateUVW(uvw, dphase, vb);
833 6036 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
834 :
835 :
836 :
837 : // Take care of translation of Bools to Integer
838 6036 : Int idopsf=0;
839 6036 : if(dopsf) idopsf=1;
840 :
841 6036 : Vector<Int> rowFlags(vb.nRows());
842 6036 : rowFlags=0;
843 6036 : rowFlags(vb.flagRow())=true;
844 6036 : if(!usezero_p) {
845 2054272 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
846 2048436 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
847 : }
848 : }
849 :
850 :
851 : Bool del;
852 : // IPosition s(flags.shape());
853 6036 : Vector<Int> s(flags.shape().nelements());
854 6036 : convertArray(s, flags.shape().asVector());
855 6036 : Int nvp=s[0];
856 6036 : Int nvc=s[1];
857 6036 : Int nvisrow=s[2];
858 :
859 6036 : Int csamp=convSampling;
860 :
861 : Bool uvwcopy;
862 6036 : const Double *uvwstor=uvw.getStorage(uvwcopy);
863 : Bool gridcopy;
864 : Bool convcopy;
865 6036 : const Complex *convstor=convFunc.getStorage(convcopy);
866 6036 : Cube<Int> loc(3, nvc, nRow);
867 6036 : Cube<Int> off(3, nvc, nRow);
868 6036 : Matrix<Complex> phasor(nvc, nRow);
869 : Bool delphase;
870 6036 : Complex * phasorstor=phasor.getStorage(delphase);
871 6036 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
872 6036 : const Double * scalestor=uvScale.getStorage(del);
873 6036 : const Double * offsetstor=uvOffset.getStorage(del);
874 : Bool delloc, deloff;
875 6036 : Int * locstor=loc.getStorage(delloc);
876 6036 : Int * offstor=off.getStorage(deloff);
877 6036 : const Double *dpstor=dphase.getStorage(del);
878 6036 : Double cinv=Double(1.0)/C::c;
879 : Int irow;
880 6036 : Int dow=1;
881 6036 : Int nth=1;
882 : #ifdef _OPENMP
883 6036 : if(numthreads_p >0){
884 0 : nth=min(numthreads_p, omp_get_max_threads());
885 : }
886 : else{
887 6036 : nth= omp_get_max_threads();
888 : }
889 : //nth=min(4,nth);
890 : #endif
891 :
892 :
893 6036 : #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 6036 : Int maxx=0;
906 6036 : Int minx=1000000;
907 6036 : Int maxy=0;
908 6036 : Int miny=1000000;
909 : //Int maxsupp=max(convSupport);
910 6036 : loc.putStorage(locstor, delloc);
911 6036 : maxx=max(loc.yzPlane(0));
912 6036 : maxx=min(nx-1, maxx-1);
913 6036 : minx=min(loc.yzPlane(0));
914 6036 : minx=max(0, minx-1);
915 6036 : maxy=max(loc.yzPlane(1));
916 6036 : maxy=min(ny-1, maxy-1);
917 6036 : miny=min(loc.yzPlane(1));
918 6036 : miny=max(0,miny-1);
919 6036 : locstor=loc.getStorage(delloc);
920 : //cerr << "dels " << delloc << " " << deloff << endl;
921 : //cerr << "LOC " << min(loc.xzPlane(0)) << " max " << loc.shape() << endl;
922 6036 : timemass_p +=tim.real();
923 : Int ixsub, iysub, icounter;
924 6036 : ixsub=1;
925 6036 : iysub=1;
926 6036 : if (nth >4){
927 3170 : ixsub=8;
928 3170 : iysub=8;
929 : }
930 : else {
931 2866 : ixsub=2;
932 2866 : iysub=2;
933 : }
934 : //nxsub=nx;
935 : //nysub=ny;
936 6036 : Int rbeg=startRow+1;
937 6036 : Int rend=endRow+1;
938 :
939 6036 : const Int* pmapstor=polMap.getStorage(del);
940 6036 : const Int *cmapstor=chanMap.getStorage(del);
941 6036 : Int nc=nchan;
942 6036 : Int np=npol;
943 : // Int nxp=nx;
944 : // Int nyp=ny;
945 6036 : minx=0;
946 6036 : miny=0;
947 6036 : Int nxp=nx;
948 6036 : Int nyp=ny;
949 6036 : Int nxcopy=nx;
950 6036 : Int nycopy=ny;
951 :
952 6036 : Int csize=convSize;
953 6036 : Int wcsize=wConvSize;
954 6036 : const Int * flagstor=flags.getStorage(del);
955 6036 : const Int * rowflagstor=rowFlags.getStorage(del);
956 6036 : const Int * suppstor=convSupport.getStorage(del);
957 6036 : tim.mark();
958 6036 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
959 220380 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
960 214344 : sumwgt[icounter].resize(sumWeight.shape());
961 214344 : sumwgt[icounter].set(0.0);
962 : }
963 6036 : if(doneThreadPartition_p < 0){
964 52 : xsect_p.resize(ixsub*iysub);
965 52 : ysect_p.resize(ixsub*iysub);
966 52 : nxsect_p.resize(ixsub*iysub);
967 52 : nysect_p.resize(ixsub*iysub);
968 1160 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
969 1108 : findGridSector(nxp, nyp, ixsub, iysub, minx, miny, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
970 : }
971 : }
972 6036 : Vector<Int> xsect, ysect, nxsect, nysect;
973 6036 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
974 :
975 6036 : if(!useDoubleGrid_p){
976 200 : Complex *gridstor=griddedData.getStorage(gridcopy);
977 200 : #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 200 : if(dopsf && (nth > 4))
1016 0 : tweakGridSector(nx, ny, ixsub, iysub);
1017 200 : timegrid_p+=tim.real();
1018 :
1019 13000 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1020 12800 : sumWeight=sumWeight+sumwgt[icounter];
1021 : }
1022 200 : griddedData.putStorage(gridstor, gridcopy);
1023 :
1024 : }
1025 : else{
1026 5836 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
1027 5836 : #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 5836 : if(dopsf && (nth > 4))
1064 918 : tweakGridSector(nx, ny, ixsub, iysub);
1065 5836 : timegrid_p+=tim.real();
1066 :
1067 207380 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1068 201544 : sumWeight=sumWeight+sumwgt[icounter];
1069 : }
1070 5836 : griddedData2.putStorage(gridstor, gridcopy);
1071 : }
1072 6036 : uvw.freeStorage(uvwstor, uvwcopy);
1073 6036 : convFunc.freeStorage(convstor, convcopy);
1074 :
1075 6036 : if(!dopsf)
1076 3883 : data.freeStorage(datStorage, isCopy);
1077 6036 : elWeight.freeStorage(wgtStorage,iswgtCopy);
1078 6036 : }
1079 :
1080 :
1081 :
1082 1730 : void WProjectFT::get(VisBuffer2& vb, Int row)
1083 : {
1084 :
1085 :
1086 1730 : findConvFunction(*image, vb);
1087 : // If row is -1 then we pass through all rows
1088 : Int startRow, endRow, nRow;
1089 1730 : if (row==-1) {
1090 1730 : nRow=vb.nRows();
1091 1730 : startRow=0;
1092 1730 : 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 1730 : matchChannel(vb);
1104 :
1105 : //No point in reading data if its not matching in frequency
1106 1730 : if(max(chanMap)==-1)
1107 0 : return;
1108 :
1109 :
1110 : // Get the uvws in a form that Fortran can use
1111 1730 : Matrix<Double> uvw(negateUV(vb));
1112 1730 : Vector<Double> dphase(vb.nRows());
1113 1730 : dphase=0.0;
1114 :
1115 1730 : rotateUVW(uvw, dphase, vb);
1116 1730 : 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 1730 : Cube<Complex> data;
1129 1730 : Cube<Int> flags;
1130 1730 : getInterpolateArrays(vb, data, flags);
1131 :
1132 :
1133 :
1134 : Complex *datStorage;
1135 : Bool isCopy;
1136 1730 : datStorage=data.getStorage(isCopy);
1137 :
1138 1730 : Vector<Int> rowFlags(vb.nRows());
1139 1730 : rowFlags=0;
1140 1730 : rowFlags(vb.flagRow())=true;
1141 1730 : if(!usezero_p) {
1142 538560 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1143 537030 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
1144 : }
1145 : }
1146 :
1147 1730 : Int nvp=data.shape()(0);
1148 1730 : Int nvc=data.shape()(1);
1149 1730 : Int nvisrow=data.shape()(2);
1150 1730 : Int nc=nchan;
1151 1730 : Int np=npol;
1152 1730 : Int nxp=nx;
1153 1730 : Int nyp=ny;
1154 1730 : Cube<Int> loc(3, nvc, nvisrow);
1155 1730 : Cube<Int> off(3, nvc, nRow);
1156 1730 : Int csamp=convSampling;
1157 1730 : Int csize=convSize;
1158 1730 : Int wcsize=wConvSize;
1159 1730 : Matrix<Complex> phasor(nvc, nRow);
1160 : Bool delphase;
1161 : Bool del;
1162 1730 : Complex * phasorstor=phasor.getStorage(delphase);
1163 1730 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1164 1730 : const Double * scalestor=uvScale.getStorage(del);
1165 1730 : const Double * offsetstor=uvOffset.getStorage(del);
1166 1730 : Int * locstor=loc.getStorage(del);
1167 1730 : Int * offstor=off.getStorage(del);
1168 1730 : const Int * flagstor=flags.getStorage(del);
1169 1730 : const Int * rowflagstor=rowFlags.getStorage(del);
1170 1730 : const Double *dpstor=dphase.getStorage(del);
1171 : Bool uvwcopy;
1172 1730 : const Double *uvwstor=uvw.getStorage(uvwcopy);
1173 : Bool gridcopy;
1174 1730 : const Complex *gridstor=griddedData.getStorage(gridcopy);
1175 : Bool convcopy;
1176 1730 : const Complex *convstor=convFunc.getStorage(convcopy);
1177 1730 : const Int* pmapstor=polMap.getStorage(del);
1178 1730 : const Int *cmapstor=chanMap.getStorage(del);
1179 1730 : const Int * suppstor=convSupport.getStorage(del);
1180 : Int irow;
1181 1730 : Int nth=1;
1182 : #ifdef _OPENMP
1183 1730 : if(numthreads_p >0){
1184 0 : nth=min(numthreads_p, omp_get_max_threads());
1185 : }
1186 : else{
1187 1730 : nth= omp_get_max_threads();
1188 : }
1189 : // nth=min(4,nth);
1190 : #endif
1191 1730 : Int dow=1;
1192 1730 : Double cinv=Double(1.0)/C::c;
1193 1730 : #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 1730 : Int rbeg=startRow+1;
1206 1730 : Int rend=endRow+1;
1207 1730 : Int npart=nth;
1208 1730 : Timer tim;
1209 1730 : tim.mark();
1210 :
1211 1730 : Int ix=0;
1212 1730 : #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 1730 : data.putStorage(datStorage, isCopy);
1244 1730 : uvw.freeStorage(uvwstor, uvwcopy);
1245 1730 : griddedData.freeStorage(gridstor, gridcopy);
1246 1730 : convFunc.freeStorage(convstor, convcopy);
1247 1730 : timedegrid_p+=tim.real();
1248 :
1249 1730 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1250 1730 : }
1251 :
1252 :
1253 :
1254 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1255 : // return the resulting image
1256 78 : ImageInterface<Complex>& WProjectFT::getImage(Matrix<Float>& weights,
1257 : Bool normalize)
1258 : {
1259 : //AlwaysAssert(lattice, AipsError);
1260 78 : AlwaysAssert(image, AipsError);
1261 :
1262 78 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
1263 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
1264 : }
1265 78 : logIO() << LogOrigin("WProjectFT", "getImage") << LogIO::NORMAL;
1266 :
1267 78 : weights.resize(sumWeight.shape());
1268 :
1269 78 : convertArray(weights, sumWeight);
1270 :
1271 : // If the weights are all zero then we cannot normalize
1272 : // otherwise we don't care.
1273 78 : 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 78 : << "Starting FFT and scaling of image" << LogIO::POST;
1288 :
1289 78 : if(useDoubleGrid_p){
1290 77 : ArrayLattice<DComplex> darrayLattice(griddedData2);
1291 77 : LatticeFFT::cfft2d(darrayLattice,false);
1292 77 : griddedData.resize(griddedData2.shape());
1293 77 : convertArray(griddedData, griddedData2);
1294 77 : SynthesisUtilMethods::getResource("mem peak in getImage");
1295 77 : griddedData2.resize();
1296 77 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1297 77 : lattice=arrayLattice;
1298 77 : }else{
1299 1 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1300 1 : lattice=arrayLattice;
1301 1 : LatticeFFT::cfft2d(*lattice,false);
1302 :
1303 : }
1304 :
1305 :
1306 78 : const IPosition latticeShape = lattice->shape();
1307 :
1308 :
1309 : {
1310 78 : Int inx = lattice->shape()(0);
1311 78 : Int iny = lattice->shape()(1);
1312 78 : Vector<Complex> correction(inx);
1313 78 : correction=Complex(1.0, 0.0);
1314 :
1315 78 : Int npixCorr= max(nx,ny);
1316 78 : Vector<Float> sincConv(npixCorr);
1317 80778 : for (Int ix=0;ix<npixCorr;ix++) {
1318 80700 : Float x=C::pi*Float(ix-npixCorr/2)/(Float(npixCorr)*Float(convSampling));
1319 80700 : if(ix==npixCorr/2) {
1320 78 : sincConv(ix)=1.0;
1321 : }
1322 : else {
1323 80622 : sincConv(ix)=sin(x)/x;
1324 : }
1325 : }
1326 :
1327 78 : IPosition cursorShape(4, inx, 1, 1, 1);
1328 78 : IPosition axisPath(4, 0, 1, 2, 3);
1329 78 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1330 78 : LatticeIterator<Complex> lix(*lattice, lsx);
1331 117298 : for(lix.reset();!lix.atEnd();lix++) {
1332 117220 : Int pol=lix.position()(2);
1333 117220 : Int chan=lix.position()(3);
1334 117220 : if(weights(pol, chan)!=0.0) {
1335 117220 : Int iy=lix.position()(1);
1336 117220 : gridder->correctX1D(correction,iy);
1337 160561820 : for (Int ix=0;ix<nx;ix++) {
1338 160444600 : correction(ix)*=sincConv(ix)*sincConv(iy);
1339 : }
1340 117220 : lix.rwVectorCursor()/=correction;
1341 117220 : if(normalize) {
1342 100 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
1343 100 : lix.rwCursor()*=rnorm;
1344 : }
1345 : else {
1346 117120 : Complex rnorm(Float(inx)*Float(iny));
1347 117120 : lix.rwCursor()*=rnorm;
1348 : }
1349 : }
1350 : else {
1351 0 : lix.woCursor()=0.0;
1352 : }
1353 : }
1354 78 : }
1355 :
1356 78 : if(!isTiled) {
1357 78 : LatticeLocker lock1 (*(image), FileLocker::Write);
1358 : // Check the section from the image BEFORE converting to a lattice
1359 156 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
1360 156 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1361 78 : IPosition stride(4, 1);
1362 156 : IPosition trc(blc+image->shape()-stride);
1363 :
1364 : // Do the copy
1365 78 : image->put(griddedData(blc, trc));
1366 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
1367 78 : griddedData.resize(IPosition(1,0));
1368 78 : }
1369 78 : }
1370 78 : if(!lattice.null()) lattice=nullptr;
1371 78 : if(!arrayLattice.null()) lattice=nullptr;
1372 78 : griddedData.resize();
1373 78 : 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 0 : }
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 0 : }
1511 :
1512 15548 : void WProjectFT::ok() {
1513 15548 : AlwaysAssert(image, AipsError);
1514 15548 : }
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 142 : String WProjectFT::name() const {
1586 :
1587 142 : 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 252 : }
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 :
|
