Line data Source code
1 : //# WProjectFT.cc: Implementation of WProjectFT class
2 : //# Copyright (C) 2003-2016
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <msvis/MSVis/VisibilityIterator2.h>
29 : #include <casacore/casa/Quanta/UnitMap.h>
30 : #include <casacore/casa/Quanta/MVTime.h>
31 : #include <casacore/casa/Quanta/UnitVal.h>
32 : #include <casacore/casa/Utilities/CountedPtr.h>
33 : #include <casacore/measures/Measures/Stokes.h>
34 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
35 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
36 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
37 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
38 : #include <casacore/coordinates/Coordinates/Projection.h>
39 : #include <casacore/ms/MeasurementSets/MSColumns.h>
40 : #include <casacore/casa/BasicSL/Constants.h>
41 : #include <casacore/scimath/Mathematics/FFTServer.h>
42 : #include <synthesis/TransformMachines2/WProjectFT.h>
43 : #include <synthesis/TransformMachines2/WPConvFunc.h>
44 : #include <casacore/scimath/Mathematics/RigidVector.h>
45 : #include <msvis/MSVis/StokesVector.h>
46 : #include <synthesis/TransformMachines/StokesImageUtil.h>
47 : #include <casacore/images/Images/ImageInterface.h>
48 : #include <casacore/images/Images/PagedImage.h>
49 : #include <casacore/casa/Containers/Block.h>
50 : #include <casacore/casa/Containers/Record.h>
51 : #include <casacore/casa/Arrays/ArrayLogical.h>
52 : #include <casacore/casa/Arrays/ArrayMath.h>
53 : #include <casacore/casa/Arrays/Array.h>
54 : #include <casacore/casa/Arrays/MaskedArray.h>
55 : #include <casacore/casa/Arrays/Vector.h>
56 : #include <casacore/casa/Arrays/Slice.h>
57 : #include <casacore/casa/Arrays/Matrix.h>
58 : #include <casacore/casa/Arrays/Cube.h>
59 : #include <casacore/casa/Arrays/MatrixIter.h>
60 : #include <casacore/casa/BasicSL/String.h>
61 : #include <casacore/casa/Utilities/Assert.h>
62 : #include <casacore/casa/Exceptions/Error.h>
63 : #include <iostream>
64 : #include <iomanip>
65 : #include <casacore/lattices/Lattices/ArrayLattice.h>
66 : #include <casacore/lattices/Lattices/SubLattice.h>
67 : #include <casacore/lattices/LRegions/LCBox.h>
68 : #include <casacore/lattices/LEL/LatticeExpr.h>
69 : #include <casacore/lattices/Lattices/LatticeCache.h>
70 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
71 : #include <casacore/lattices/Lattices/LatticeIterator.h>
72 : #include <casacore/lattices/Lattices/LatticeStepper.h>
73 : #include <casacore/casa/Utilities/CompositeNumber.h>
74 : #include <casacore/casa/OS/Timer.h>
75 : #include <casacore/casa/OS/HostInfo.h>
76 : #include <sstream>
77 : #ifdef _OPENMP
78 : #include <omp.h>
79 : #endif
80 : using namespace casacore;
81 : namespace casa { //# NAMESPACE CASA - BEGIN
82 : namespace refim { //#namespace for imaging refactoring
83 :
84 : using namespace casacore;
85 : using namespace casa;
86 : using namespace casacore;
87 : using namespace casa::vi;
88 : using namespace casacore;
89 : using namespace casa::refim;
90 :
91 0 : WProjectFT::WProjectFT( Int nWPlanes, Long icachesize, Int itilesize,
92 0 : Bool usezero, Bool useDoublePrec, const Double minW, const Double maxW, const Double rmsW)
93 : : FTMachine(), padding_p(1.0), nWPlanes_p(nWPlanes),
94 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
95 : gridder(0), isTiled(false),
96 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)), usezero_p(usezero),
97 0 : machineName_p("WProjectFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), minW_p(minW), maxW_p(maxW), rmsW_p(rmsW)
98 : {
99 0 : convSize=0;
100 0 : tangentSpecified_p=false;
101 0 : lastIndex_p=0;
102 0 : useDoubleGrid_p=useDoublePrec;
103 0 : wpConvFunc_p=new WPConvFunc(minW, maxW, rmsW);
104 :
105 0 : }
106 0 : WProjectFT::WProjectFT(Int nWPlanes,
107 : MPosition mLocation,
108 : Long icachesize, Int itilesize,
109 0 : Bool usezero, Float padding, Bool useDoublePrec, const Double minW, const Double maxW, const Double rmsW)
110 : : FTMachine(), padding_p(padding), nWPlanes_p(nWPlanes),
111 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
112 : gridder(0), isTiled(false),
113 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
114 : usezero_p(usezero),
115 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)
116 :
117 : {
118 0 : convSize=0;
119 0 : savedWScale_p=0.0;
120 0 : tangentSpecified_p=false;
121 0 : mLocation_p=mLocation;
122 0 : lastIndex_p=0;
123 0 : wpConvFunc_p=new WPConvFunc(minW, maxW, rmsW);
124 0 : useDoubleGrid_p=useDoublePrec;
125 0 : }
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 : : FTMachine(), padding_p(padding), nWPlanes_p(nWPlanes),
132 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
133 : gridder(0), isTiled(false),
134 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
135 : usezero_p(usezero),
136 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 0 : 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 0 : isTiled=false;
226 0 : CompositeNumber cn(uInt(image->shape()(0)*2));
227 0 : nx = cn.nextLargerEven(Int(padding_p*Float(image->shape()(0))-0.5));
228 0 : ny = cn.nextLargerEven(Int(padding_p*Float(image->shape()(1))-0.5));
229 0 : npol = image->shape()(2);
230 0 : 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 0 : isTiled=false;
241 :
242 :
243 0 : sumWeight.resize(npol, nchan);
244 :
245 0 : wConvSize=max(0, nWPlanes_p);
246 0 : convSupport.resize(wConvSize);
247 0 : convSupport=0;
248 :
249 0 : uvScale.resize(3);
250 0 : uvScale=0.0;
251 0 : uvScale(0)=Float(nx)*image->coordinates().increment()(0);
252 0 : uvScale(1)=Float(ny)*image->coordinates().increment()(1);
253 0 : if(savedWScale_p==0.0){
254 0 : uvScale(2)=Float(wConvSize)*abs(image->coordinates().increment()(0));
255 : }
256 : else{
257 0 : uvScale(2)=savedWScale_p;
258 : }
259 0 : uvOffset.resize(3);
260 0 : uvOffset(0)=nx/2;
261 0 : uvOffset(1)=ny/2;
262 0 : uvOffset(2)=0;
263 :
264 :
265 :
266 0 : if(gridder) delete gridder; gridder=0;
267 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
268 0 : uvScale, uvOffset,
269 0 : "SF");
270 :
271 : // Set up image cache needed for gridding.
272 0 : 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 0 : if(isTiled) {
277 0 : Float tileOverlap=0.5;
278 0 : tilesize=min(256,tilesize);
279 0 : IPosition tileShape=IPosition(4,tilesize,tilesize,npol,nchan);
280 0 : Vector<Float> tileOverlapVec(4);
281 0 : tileOverlapVec=0.0;
282 0 : tileOverlapVec(0)=tileOverlap;
283 0 : tileOverlapVec(1)=tileOverlap;
284 0 : Int tmpCacheVal=static_cast<Int>(cachesize);
285 0 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
286 : tileOverlapVec,
287 0 : (tileOverlap>0.0));
288 :
289 : }
290 0 : }
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 0 : void WProjectFT::findConvFunction(const ImageInterface<Complex>& image,
322 : const VisBuffer2& vb) {
323 :
324 :
325 :
326 :
327 :
328 0 : wpConvFunc_p->findConvFunction(image, vb, wConvSize, uvScale, uvOffset,
329 0 : padding_p,
330 0 : convSampling,
331 0 : convFunc, convSize, convSupport,
332 0 : savedWScale_p);
333 :
334 0 : nWPlanes_p=convSupport.nelements();
335 0 : wConvSize=max(1,nWPlanes_p);
336 0 : uvScale(2)=savedWScale_p;
337 :
338 :
339 :
340 0 : }
341 :
342 0 : void WProjectFT::initializeToVis(ImageInterface<Complex>& iimage,
343 : const VisBuffer2& vb)
344 : {
345 0 : image=&iimage;
346 0 : toVis_p=true;
347 0 : ok();
348 :
349 : // if(convSize==0) {
350 0 : init();
351 : // }
352 0 : findConvFunction(*image, vb);
353 :
354 :
355 : // Initialize the maps for polarization and channel. These maps
356 : // translate visibility indices into image indices
357 0 : 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 0 : 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 0 : prepGridForDegrid();
377 0 : }
378 :
379 0 : void WProjectFT::prepGridForDegrid(){
380 0 : IPosition gridShape(4, nx, ny, npol, nchan);
381 0 : griddedData.resize(gridShape);
382 0 : griddedData=Complex(0.0);
383 :
384 :
385 0 : IPosition stride(4, 1);
386 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
387 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
388 0 : IPosition trc(blc+image->shape()-stride);
389 :
390 0 : IPosition start(4, 0);
391 0 : griddedData(blc, trc) = image->getSlice(start, image->shape());
392 :
393 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
394 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
395 0 : lattice=arrayLattice;
396 :
397 0 : logIO() << LogIO::DEBUGGING << "Starting FFT of image" << LogIO::POST;
398 :
399 :
400 0 : Vector<Float> sincConvX(nx);
401 0 : for (Int ix=0;ix<nx;ix++) {
402 0 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
403 0 : if(ix==nx/2) {
404 0 : sincConvX(ix)=1.0;
405 : }
406 : else {
407 0 : sincConvX(ix)=sin(x)/x;
408 : }
409 : }
410 0 : Vector<Float> sincConvY(ny);
411 0 : for (Int ix=0;ix<ny;ix++) {
412 0 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
413 0 : if(ix==ny/2) {
414 0 : sincConvY(ix)=1.0;
415 : }
416 : else {
417 0 : sincConvY(ix)=sin(x)/x;
418 : }
419 : }
420 :
421 0 : Vector<Complex> correction(nx);
422 0 : correction=Complex(1.0, 0.0);
423 : // Do the Grid-correction
424 0 : IPosition cursorShape(4, nx, 1, 1, 1);
425 0 : IPosition axisPath(4, 0, 1, 2, 3);
426 0 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
427 0 : LatticeIterator<Complex> lix(*lattice, lsx);
428 0 : for(lix.reset();!lix.atEnd();lix++) {
429 0 : Int iy=lix.position()(1);
430 0 : gridder->correctX1D(correction,iy);
431 0 : for (Int ix=0;ix<nx;ix++) {
432 0 : correction(ix)/=(sincConvX(ix)*sincConvY(iy));
433 : }
434 0 : lix.rwVectorCursor()/=correction;
435 : }
436 :
437 : // Now do the FFT2D in place
438 0 : LatticeFFT::cfft2d(*lattice);
439 :
440 0 : logIO() << LogIO::DEBUGGING << "Finished FFT" << LogIO::POST;
441 :
442 :
443 :
444 0 : }
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 : }
466 0 : }
467 :
468 :
469 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
470 : // grid.
471 0 : void WProjectFT::initializeToSky(ImageInterface<Complex>& iimage,
472 : Matrix<Float>& weight,
473 : const VisBuffer2& vb)
474 : {
475 : // image always points to the image
476 0 : image=&iimage;
477 0 : toVis_p=false;
478 :
479 : // if(convSize==0) {
480 0 : init();
481 : // }
482 0 : findConvFunction(*image, vb);
483 :
484 :
485 : // Initialize the maps for polarization and channel. These maps
486 : // translate visibility indices into image indices
487 0 : 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 0 : isTiled=false;
501 0 : sumWeight=0.0;
502 0 : weight.resize(sumWeight.shape());
503 0 : 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 0 : 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 0 : IPosition gridShape(4, nx, ny, npol, nchan);
515 0 : if(!useDoubleGrid_p){
516 0 : griddedData.resize(gridShape);
517 0 : 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 : }
527 : //AlwaysAssert(lattice, AipsError);
528 :
529 0 : }
530 :
531 0 : void WProjectFT::finalizeToSky()
532 : {
533 0 : logIO() << LogOrigin("WProjectFT", "finalizeToSky") << LogIO::NORMAL;
534 0 : logIO() <<LogIO::NORMAL2<< "Time to massage data " << timemass_p << LogIO::POST;
535 0 : logIO() << LogIO::NORMAL2 <<"Time to grid data " << timegrid_p << LogIO::POST;
536 0 : timemass_p=0.0;
537 0 : 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 0 : if(isTiled) {
541 0 : logIO() << LogOrigin("WProjectFT", "finalizeToSky") << LogIO::NORMAL;
542 :
543 0 : AlwaysAssert(image, AipsError);
544 0 : AlwaysAssert(imageCache, AipsError);
545 0 : imageCache->flush();
546 0 : ostringstream o;
547 0 : imageCache->showCacheStatistics(o);
548 0 : logIO() << o.str() << LogIO::POST;
549 : }
550 0 : }
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 0 : 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 0 : 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 0 : if(max(chanMap)==-1)
782 0 : return;
783 0 : Timer tim;
784 0 : tim.mark();
785 :
786 : //const Matrix<Float> *imagingweight;
787 : //imagingweight=&(vb.imagingWeight());
788 0 : Matrix<Float> imagingweight;
789 0 : getImagingWeight(imagingweight, vb);
790 0 : if(dopsf) type=FTMachine::PSF;
791 :
792 0 : Cube<Complex> data;
793 : //Fortran gridder need the flag as ints
794 0 : Cube<Int> flags;
795 0 : Matrix<Float> elWeight;
796 0 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
797 :
798 :
799 : Bool iswgtCopy;
800 : const Float *wgtStorage;
801 0 : wgtStorage=elWeight.getStorage(iswgtCopy);
802 :
803 :
804 : Bool isCopy;
805 0 : const Complex *datStorage=0;
806 :
807 0 : if(!dopsf)
808 0 : datStorage=data.getStorage(isCopy);
809 :
810 :
811 : // If row is -1 then we pass through all rows
812 : Int startRow, endRow, nRow;
813 0 : if (row==-1) {
814 0 : nRow=vb.nRows();
815 0 : startRow=0;
816 0 : 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 0 : Matrix<Double> uvw(negateUV(vb));
829 0 : Vector<Double> dphase(vb.nRows());
830 0 : dphase=0.0;
831 :
832 0 : rotateUVW(uvw, dphase, vb);
833 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
834 :
835 :
836 :
837 : // Take care of translation of Bools to Integer
838 0 : Int idopsf=0;
839 0 : if(dopsf) idopsf=1;
840 :
841 0 : Vector<Int> rowFlags(vb.nRows());
842 0 : rowFlags=0;
843 0 : rowFlags(vb.flagRow())=true;
844 0 : 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 0 : Vector<Int> s(flags.shape().nelements());
854 0 : convertArray(s, flags.shape().asVector());
855 0 : Int nvp=s[0];
856 0 : Int nvc=s[1];
857 0 : Int nvisrow=s[2];
858 :
859 0 : Int csamp=convSampling;
860 :
861 : Bool uvwcopy;
862 0 : const Double *uvwstor=uvw.getStorage(uvwcopy);
863 : Bool gridcopy;
864 : Bool convcopy;
865 0 : const Complex *convstor=convFunc.getStorage(convcopy);
866 0 : Cube<Int> loc(3, nvc, nRow);
867 0 : Cube<Int> off(3, nvc, nRow);
868 0 : Matrix<Complex> phasor(nvc, nRow);
869 : Bool delphase;
870 0 : Complex * phasorstor=phasor.getStorage(delphase);
871 0 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
872 0 : const Double * scalestor=uvScale.getStorage(del);
873 0 : const Double * offsetstor=uvOffset.getStorage(del);
874 : Bool delloc, deloff;
875 0 : Int * locstor=loc.getStorage(delloc);
876 0 : Int * offstor=off.getStorage(deloff);
877 0 : const Double *dpstor=dphase.getStorage(del);
878 0 : Double cinv=Double(1.0)/C::c;
879 : Int irow;
880 0 : Int dow=1;
881 0 : Int nth=1;
882 : #ifdef _OPENMP
883 0 : if(numthreads_p >0){
884 0 : nth=min(numthreads_p, omp_get_max_threads());
885 : }
886 : else{
887 0 : nth= omp_get_max_threads();
888 : }
889 : //nth=min(4,nth);
890 : #endif
891 :
892 :
893 0 : #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 0 : Int maxx=0;
906 0 : Int minx=1000000;
907 0 : Int maxy=0;
908 0 : Int miny=1000000;
909 : //Int maxsupp=max(convSupport);
910 0 : loc.putStorage(locstor, delloc);
911 0 : maxx=max(loc.yzPlane(0));
912 0 : maxx=min(nx-1, maxx-1);
913 0 : minx=min(loc.yzPlane(0));
914 0 : minx=max(0, minx-1);
915 0 : maxy=max(loc.yzPlane(1));
916 0 : maxy=min(ny-1, maxy-1);
917 0 : miny=min(loc.yzPlane(1));
918 0 : miny=max(0,miny-1);
919 0 : locstor=loc.getStorage(delloc);
920 : //cerr << "dels " << delloc << " " << deloff << endl;
921 : //cerr << "LOC " << min(loc.xzPlane(0)) << " max " << loc.shape() << endl;
922 0 : timemass_p +=tim.real();
923 : Int ixsub, iysub, icounter;
924 0 : ixsub=1;
925 0 : iysub=1;
926 0 : if (nth >4){
927 0 : ixsub=8;
928 0 : iysub=8;
929 : }
930 : else {
931 0 : ixsub=2;
932 0 : iysub=2;
933 : }
934 : //nxsub=nx;
935 : //nysub=ny;
936 0 : Int rbeg=startRow+1;
937 0 : Int rend=endRow+1;
938 :
939 0 : const Int* pmapstor=polMap.getStorage(del);
940 0 : const Int *cmapstor=chanMap.getStorage(del);
941 0 : Int nc=nchan;
942 0 : Int np=npol;
943 : // Int nxp=nx;
944 : // Int nyp=ny;
945 0 : minx=0;
946 0 : miny=0;
947 0 : Int nxp=nx;
948 0 : Int nyp=ny;
949 0 : Int nxcopy=nx;
950 0 : Int nycopy=ny;
951 :
952 0 : Int csize=convSize;
953 0 : Int wcsize=wConvSize;
954 0 : const Int * flagstor=flags.getStorage(del);
955 0 : const Int * rowflagstor=rowFlags.getStorage(del);
956 0 : const Int * suppstor=convSupport.getStorage(del);
957 0 : tim.mark();
958 0 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
959 0 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
960 0 : sumwgt[icounter].resize(sumWeight.shape());
961 0 : sumwgt[icounter].set(0.0);
962 : }
963 0 : if(doneThreadPartition_p < 0){
964 0 : xsect_p.resize(ixsub*iysub);
965 0 : ysect_p.resize(ixsub*iysub);
966 0 : nxsect_p.resize(ixsub*iysub);
967 0 : nysect_p.resize(ixsub*iysub);
968 0 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
969 0 : findGridSector(nxp, nyp, ixsub, iysub, minx, miny, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
970 : }
971 : }
972 0 : Vector<Int> xsect, ysect, nxsect, nysect;
973 0 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
974 :
975 0 : 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 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 0 : uvw.freeStorage(uvwstor, uvwcopy);
1073 0 : convFunc.freeStorage(convstor, convcopy);
1074 :
1075 0 : if(!dopsf)
1076 0 : data.freeStorage(datStorage, isCopy);
1077 0 : elWeight.freeStorage(wgtStorage,iswgtCopy);
1078 : }
1079 :
1080 :
1081 :
1082 0 : void WProjectFT::get(VisBuffer2& vb, Int row)
1083 : {
1084 :
1085 :
1086 0 : findConvFunction(*image, vb);
1087 : // If row is -1 then we pass through all rows
1088 : Int startRow, endRow, nRow;
1089 0 : if (row==-1) {
1090 0 : nRow=vb.nRows();
1091 0 : startRow=0;
1092 0 : 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 0 : matchChannel(vb);
1104 :
1105 : //No point in reading data if its not matching in frequency
1106 0 : if(max(chanMap)==-1)
1107 0 : return;
1108 :
1109 :
1110 : // Get the uvws in a form that Fortran can use
1111 0 : Matrix<Double> uvw(negateUV(vb));
1112 0 : Vector<Double> dphase(vb.nRows());
1113 0 : dphase=0.0;
1114 :
1115 0 : rotateUVW(uvw, dphase, vb);
1116 0 : 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 0 : Cube<Complex> data;
1129 0 : Cube<Int> flags;
1130 0 : getInterpolateArrays(vb, data, flags);
1131 :
1132 :
1133 :
1134 : Complex *datStorage;
1135 : Bool isCopy;
1136 0 : datStorage=data.getStorage(isCopy);
1137 :
1138 0 : Vector<Int> rowFlags(vb.nRows());
1139 0 : rowFlags=0;
1140 0 : rowFlags(vb.flagRow())=true;
1141 0 : 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 0 : Int nvp=data.shape()(0);
1148 0 : Int nvc=data.shape()(1);
1149 0 : Int nvisrow=data.shape()(2);
1150 0 : Int nc=nchan;
1151 0 : Int np=npol;
1152 0 : Int nxp=nx;
1153 0 : Int nyp=ny;
1154 0 : Cube<Int> loc(3, nvc, nvisrow);
1155 0 : Cube<Int> off(3, nvc, nRow);
1156 0 : Int csamp=convSampling;
1157 0 : Int csize=convSize;
1158 0 : Int wcsize=wConvSize;
1159 0 : Matrix<Complex> phasor(nvc, nRow);
1160 : Bool delphase;
1161 : Bool del;
1162 0 : Complex * phasorstor=phasor.getStorage(delphase);
1163 0 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1164 0 : const Double * scalestor=uvScale.getStorage(del);
1165 0 : const Double * offsetstor=uvOffset.getStorage(del);
1166 0 : Int * locstor=loc.getStorage(del);
1167 0 : Int * offstor=off.getStorage(del);
1168 0 : const Int * flagstor=flags.getStorage(del);
1169 0 : const Int * rowflagstor=rowFlags.getStorage(del);
1170 0 : const Double *dpstor=dphase.getStorage(del);
1171 : Bool uvwcopy;
1172 0 : const Double *uvwstor=uvw.getStorage(uvwcopy);
1173 : Bool gridcopy;
1174 0 : const Complex *gridstor=griddedData.getStorage(gridcopy);
1175 : Bool convcopy;
1176 0 : const Complex *convstor=convFunc.getStorage(convcopy);
1177 0 : const Int* pmapstor=polMap.getStorage(del);
1178 0 : const Int *cmapstor=chanMap.getStorage(del);
1179 0 : const Int * suppstor=convSupport.getStorage(del);
1180 : Int irow;
1181 0 : Int nth=1;
1182 : #ifdef _OPENMP
1183 0 : if(numthreads_p >0){
1184 0 : nth=min(numthreads_p, omp_get_max_threads());
1185 : }
1186 : else{
1187 0 : nth= omp_get_max_threads();
1188 : }
1189 : // nth=min(4,nth);
1190 : #endif
1191 0 : Int dow=1;
1192 0 : Double cinv=Double(1.0)/C::c;
1193 0 : #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 0 : Int rbeg=startRow+1;
1206 0 : Int rend=endRow+1;
1207 0 : Int npart=nth;
1208 0 : Timer tim;
1209 0 : tim.mark();
1210 :
1211 0 : Int ix=0;
1212 0 : #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 0 : data.putStorage(datStorage, isCopy);
1244 0 : uvw.freeStorage(uvwstor, uvwcopy);
1245 0 : griddedData.freeStorage(gridstor, gridcopy);
1246 0 : convFunc.freeStorage(convstor, convcopy);
1247 0 : timedegrid_p+=tim.real();
1248 :
1249 0 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1250 : }
1251 :
1252 :
1253 :
1254 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1255 : // return the resulting image
1256 0 : ImageInterface<Complex>& WProjectFT::getImage(Matrix<Float>& weights,
1257 : Bool normalize)
1258 : {
1259 : //AlwaysAssert(lattice, AipsError);
1260 0 : AlwaysAssert(image, AipsError);
1261 :
1262 0 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
1263 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
1264 : }
1265 0 : logIO() << LogOrigin("WProjectFT", "getImage") << LogIO::NORMAL;
1266 :
1267 0 : weights.resize(sumWeight.shape());
1268 :
1269 0 : convertArray(weights, sumWeight);
1270 :
1271 : // If the weights are all zero then we cannot normalize
1272 : // otherwise we don't care.
1273 0 : 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 0 : << "Starting FFT and scaling of image" << LogIO::POST;
1288 :
1289 0 : 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 : }else{
1299 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1300 0 : lattice=arrayLattice;
1301 0 : LatticeFFT::cfft2d(*lattice,false);
1302 :
1303 : }
1304 :
1305 :
1306 0 : const IPosition latticeShape = lattice->shape();
1307 :
1308 :
1309 : {
1310 0 : Int inx = lattice->shape()(0);
1311 0 : Int iny = lattice->shape()(1);
1312 0 : Vector<Complex> correction(inx);
1313 0 : correction=Complex(1.0, 0.0);
1314 :
1315 0 : Int npixCorr= max(nx,ny);
1316 0 : Vector<Float> sincConv(npixCorr);
1317 0 : for (Int ix=0;ix<npixCorr;ix++) {
1318 0 : Float x=C::pi*Float(ix-npixCorr/2)/(Float(npixCorr)*Float(convSampling));
1319 0 : if(ix==npixCorr/2) {
1320 0 : sincConv(ix)=1.0;
1321 : }
1322 : else {
1323 0 : sincConv(ix)=sin(x)/x;
1324 : }
1325 : }
1326 :
1327 0 : IPosition cursorShape(4, inx, 1, 1, 1);
1328 0 : IPosition axisPath(4, 0, 1, 2, 3);
1329 0 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1330 0 : LatticeIterator<Complex> lix(*lattice, lsx);
1331 0 : for(lix.reset();!lix.atEnd();lix++) {
1332 0 : Int pol=lix.position()(2);
1333 0 : Int chan=lix.position()(3);
1334 0 : if(weights(pol, chan)!=0.0) {
1335 0 : Int iy=lix.position()(1);
1336 0 : gridder->correctX1D(correction,iy);
1337 0 : for (Int ix=0;ix<nx;ix++) {
1338 0 : correction(ix)*=sincConv(ix)*sincConv(iy);
1339 : }
1340 0 : lix.rwVectorCursor()/=correction;
1341 0 : if(normalize) {
1342 0 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
1343 0 : 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 : }
1355 :
1356 0 : if(!isTiled) {
1357 0 : LatticeLocker lock1 (*(image), FileLocker::Write);
1358 : // Check the section from the image BEFORE converting to a lattice
1359 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
1360 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1361 0 : IPosition stride(4, 1);
1362 0 : IPosition trc(blc+image->shape()-stride);
1363 :
1364 : // Do the copy
1365 0 : image->put(griddedData(blc, trc));
1366 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
1367 0 : griddedData.resize(IPosition(1,0));
1368 : }
1369 : }
1370 0 : if(!lattice.null()) lattice=nullptr;
1371 0 : if(!arrayLattice.null()) lattice=nullptr;
1372 0 : griddedData.resize();
1373 0 : return *image;
1374 : }
1375 :
1376 : // Get weight image
1377 0 : void WProjectFT::getWeightImage(ImageInterface<Float>& weightImage,
1378 : Matrix<Float>& weights)
1379 : {
1380 :
1381 0 : logIO() << LogOrigin("WProjectFT", "getWeightImage") << LogIO::NORMAL;
1382 :
1383 0 : weights.resize(sumWeight.shape());
1384 0 : convertArray(weights,sumWeight);
1385 :
1386 0 : const IPosition latticeShape = weightImage.shape();
1387 :
1388 0 : Int inx=latticeShape(0);
1389 0 : Int iny=latticeShape(1);
1390 :
1391 0 : IPosition loc(2, 0);
1392 0 : IPosition cursorShape(4, inx, iny, 1, 1);
1393 0 : IPosition axisPath(4, 0, 1, 2, 3);
1394 0 : LatticeStepper lsx(latticeShape, cursorShape, axisPath);
1395 0 : LatticeIterator<Float> lix(weightImage, lsx);
1396 0 : for(lix.reset();!lix.atEnd();lix++) {
1397 0 : Int pol=lix.position()(2);
1398 0 : Int chan=lix.position()(3);
1399 0 : lix.rwCursor()=weights(pol,chan);
1400 : }
1401 0 : }
1402 :
1403 0 : Bool WProjectFT::toRecord(String& error,
1404 : RecordInterface& outRec, Bool withImage, const String diskimage)
1405 : {
1406 :
1407 : /*
1408 :
1409 : CountedPtr<WPConvFunc> wpConvFunc_p;
1410 : */
1411 :
1412 : // Save the current WProjectFT object to an output state record
1413 0 : Bool retval = true;
1414 : //save the base class variables
1415 : //this is a memory hog and slow on saving and recovering...better to recompute convfunctions
1416 : /* Record wpconvrec;
1417 : if(wpConvFunc_p->toRecord(wpconvrec))
1418 : outRec.defineRecord("wpconvfunc", wpconvrec);
1419 : */
1420 0 : Float elpadd=padding_p;
1421 0 : if(toVis_p && withImage)
1422 0 : elpadd=1.0;
1423 0 : if(!FTMachine::toRecord(error, outRec, withImage, diskimage))
1424 0 : return false;
1425 :
1426 0 : outRec.define("uvscale", uvScale);
1427 0 : outRec.define("uvoffset", uvOffset);
1428 0 : outRec.define("istiled", isTiled);
1429 0 : outRec.define("cachesize", Int64(cachesize));
1430 0 : outRec.define("tilesize", tilesize);
1431 0 : outRec.define("maxabsdata", maxAbsData);
1432 0 : Vector<Int> center_loc(4), offset_loc(4);
1433 0 : for (Int k=0; k<4 ; k++){
1434 0 : center_loc(k)=centerLoc(k);
1435 0 : offset_loc(k)=offsetLoc(k);
1436 : }
1437 0 : outRec.define("centerloc", center_loc);
1438 0 : outRec.define("offsetloc", offset_loc);
1439 0 : outRec.define("padding", elpadd);
1440 0 : outRec.define("nwplanes", nWPlanes_p);
1441 0 : outRec.define("savedwscale", savedWScale_p);
1442 0 : outRec.define("usezero", usezero_p);
1443 : ///no need to save convfunc as it can be big and is recalculated anyways
1444 : ///outRec.define("convfunc", convFunc);
1445 0 : outRec.define("convsampling", convSampling);
1446 0 : outRec.define("convsize", convSize);
1447 0 : outRec.define("convsupport", convSupport);
1448 0 : outRec.define("convsizes", convSizes_p);
1449 0 : outRec.define("wconvsize", wConvSize);
1450 0 : outRec.define("lastindex", lastIndex_p);
1451 0 : outRec.define("minw", minW_p);
1452 0 : outRec.define("maxw", maxW_p);
1453 0 : outRec.define("rmsw", rmsW_p);
1454 :
1455 :
1456 :
1457 0 : return retval;
1458 : }
1459 :
1460 0 : Bool WProjectFT::fromRecord(String& error,
1461 : const RecordInterface& inRec)
1462 : {
1463 0 : if(!FTMachine::fromRecord(error, inRec))
1464 0 : return false;
1465 0 : machineName_p="WProjectFT";
1466 0 : Bool retval = true;
1467 0 : imageCache=0; lattice=0; arrayLattice=0;
1468 0 : inRec.get("uvscale", uvScale);
1469 0 : inRec.get("uvoffset", uvOffset);
1470 0 : inRec.get("istiled", isTiled);
1471 0 : cachesize=inRec.asInt64("cachesize");
1472 0 : inRec.get("tilesize", tilesize);
1473 0 : inRec.get("maxabsdata", maxAbsData);
1474 0 : Vector<Int> center_loc(4), offset_loc(4);
1475 0 : inRec.get("centerloc", center_loc);
1476 0 : inRec.get("offsetloc", offset_loc);
1477 0 : uInt ndim4 = 4;
1478 0 : centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2),
1479 0 : center_loc(3));
1480 0 : offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2),
1481 0 : offset_loc(3));
1482 0 : inRec.get("minw", minW_p);
1483 0 : inRec.get("maxw", maxW_p);
1484 0 : inRec.get("rmsw", rmsW_p);
1485 0 : if(inRec.isDefined("wpconvfunc")){
1486 0 : wpConvFunc_p=new WPConvFunc(inRec.asRecord("wpconvfunc"));
1487 : }
1488 : else{
1489 0 : wpConvFunc_p=new WPConvFunc(minW_p, maxW_p, rmsW_p);
1490 : }
1491 0 : inRec.get("padding", padding_p);
1492 0 : inRec.get("nwplanes", nWPlanes_p);
1493 0 : inRec.get("savedwscale", savedWScale_p);
1494 0 : inRec.get("usezero", usezero_p);
1495 : //inRec.get("convfunc", convFunc);
1496 0 : convFunc.resize();
1497 0 : inRec.get("convsampling", convSampling);
1498 0 : inRec.get("convsize", convSize);
1499 0 : inRec.get("convsupport", convSupport);
1500 0 : inRec.get("convsizes", convSizes_p);
1501 0 : inRec.get("wconvsize", wConvSize);
1502 0 : inRec.get("lastindex", lastIndex_p);
1503 0 : gridder=0;
1504 : ///setup some of the parameters
1505 0 : init();
1506 :
1507 :
1508 :
1509 0 : return retval;
1510 : }
1511 :
1512 0 : void WProjectFT::ok() {
1513 0 : AlwaysAssert(image, AipsError);
1514 0 : }
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 : }
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 :
|
