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