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