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