Line data Source code
1 : //# GridFT.cc: Implementation of GridFT class
2 : //# Copyright (C) 1997-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 Library 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/UnitVal.h>
31 : #include <casacore/measures/Measures/Stokes.h>
32 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
33 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
34 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
35 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
36 : #include <casacore/coordinates/Coordinates/Projection.h>
37 : #include <casacore/ms/MeasurementSets/MSColumns.h>
38 : #include <casacore/casa/BasicSL/Constants.h>
39 : #include <casacore/scimath/Mathematics/FFTServer.h>
40 : #include <casacore/scimath/Mathematics/RigidVector.h>
41 : #include <msvis/MSVis/StokesVector.h>
42 : #include <synthesis/TransformMachines/StokesImageUtil.h>
43 : #include <synthesis/TransformMachines2/GridFT.h>
44 : #include <synthesis/Utilities/FFT2D.h>
45 : #include <msvis/MSVis/VisBuffer2.h>
46 : #include <casacore/images/Images/ImageInterface.h>
47 : #include <casacore/images/Images/PagedImage.h>
48 : #include <casacore/casa/Containers/Block.h>
49 : #include <casacore/casa/Containers/Record.h>
50 : #include <casacore/casa/Arrays/ArrayLogical.h>
51 : #include <casacore/casa/Arrays/ArrayMath.h>
52 : #include <casacore/casa/Arrays/Array.h>
53 : #include <casacore/casa/Arrays/MaskedArray.h>
54 : #include <casacore/casa/Arrays/Vector.h>
55 : #include <casacore/casa/Arrays/Matrix.h>
56 : #include <casacore/casa/Arrays/Cube.h>
57 : #include <casacore/casa/Arrays/MatrixIter.h>
58 : #include <casacore/casa/BasicSL/String.h>
59 : #include <casacore/casa/Utilities/Assert.h>
60 : #include <casacore/casa/Exceptions/Error.h>
61 : #include <casacore/lattices/Lattices/ArrayLattice.h>
62 : #include <casacore/measures/Measures/UVWMachine.h>
63 : #include <casacore/lattices/Lattices/SubLattice.h>
64 : #include <casacore/lattices/LRegions/LCBox.h>
65 : #include <casacore/lattices/Lattices/LatticeCache.h>
66 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
67 : #include <casacore/lattices/Lattices/LatticeIterator.h>
68 : #include <casacore/lattices/Lattices/LatticeStepper.h>
69 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
70 : #include <casacore/casa/Utilities/CompositeNumber.h>
71 : #include <casacore/casa/OS/Timer.h>
72 : #include <sstream>
73 : #ifdef _OPENMP
74 : #include <omp.h>
75 : #endif
76 :
77 : using namespace casacore;
78 : namespace casa { //# NAMESPACE CASA - BEGIN
79 : namespace refim {//# namespace for imaging refactor
80 : using namespace casacore;
81 : using namespace casa;
82 : using namespace casacore;
83 : using namespace casa::refim;
84 :
85 0 : GridFT::GridFT() : FTMachine(), padding_p(1.0), imageCache(0), cachesize(1000000), tilesize(1000), gridder(0), isTiled(false), convType("SF"),
86 0 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
87 0 : usezero_p(false), noPadding_p(false), usePut2_p(false),
88 0 : machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0){
89 :
90 0 : }
91 21 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType, Float padding,
92 21 : Bool usezero, Bool useDoublePrec)
93 21 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize), tilesize(itilesize),
94 21 : gridder(0), isTiled(false), convType(iconvType),
95 21 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
96 21 : usezero_p(usezero), noPadding_p(false), usePut2_p(false),
97 42 : machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
98 : {
99 21 : useDoubleGrid_p=useDoublePrec;
100 : // peek=NULL;
101 21 : }
102 :
103 2479 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
104 2479 : MPosition mLocation, Float padding, Bool usezero, Bool useDoublePrec)
105 2479 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
106 2479 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
107 2479 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
108 4958 : usePut2_p(false), machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
109 : {
110 2479 : mLocation_p=mLocation;
111 2479 : tangentSpecified_p=false;
112 2479 : useDoubleGrid_p=useDoublePrec;
113 : // peek=NULL;
114 2479 : }
115 :
116 0 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
117 0 : MDirection mTangent, Float padding, Bool usezero, Bool useDoublePrec)
118 0 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
119 0 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
120 0 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
121 0 : usePut2_p(false), machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
122 : {
123 0 : mTangent_p=mTangent;
124 0 : tangentSpecified_p=true;
125 0 : useDoubleGrid_p=useDoublePrec;
126 : // peek=NULL;
127 0 : }
128 :
129 7 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
130 : MPosition mLocation, MDirection mTangent, Float padding,
131 7 : Bool usezero, Bool useDoublePrec)
132 7 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
133 7 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
134 7 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
135 14 : usePut2_p(false),machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
136 : {
137 7 : mLocation_p=mLocation;
138 7 : mTangent_p=mTangent;
139 7 : tangentSpecified_p=true;
140 7 : useDoubleGrid_p=useDoublePrec;
141 : // peek=NULL;
142 7 : }
143 :
144 0 : GridFT::GridFT(const RecordInterface& stateRec)
145 0 : : FTMachine()
146 : {
147 : // Construct from the input state record
148 0 : String error;
149 0 : if (!fromRecord(error, stateRec)) {
150 0 : throw (AipsError("Failed to create gridder: " + error));
151 : };
152 0 : timemass_p=0.0;
153 0 : timegrid_p=0.0;
154 0 : timedegrid_p=0.0;
155 : // peek=NULL;
156 0 : }
157 :
158 : //----------------------------------------------------------------------
159 264 : GridFT& GridFT::operator=(const GridFT& other)
160 : {
161 264 : if(this!=&other) {
162 : //Do the base parameters
163 264 : FTMachine::operator=(other);
164 :
165 : //private params
166 264 : imageCache=other.imageCache;
167 264 : cachesize=other.cachesize;
168 264 : tilesize=other.tilesize;
169 264 : convType=other.convType;
170 264 : convType.upcase();
171 264 : uvScale.resize();
172 264 : uvOffset.resize();
173 264 : uvScale.assign(other.uvScale);
174 264 : uvOffset.assign(other.uvOffset);
175 264 : if(other.gridder==0)
176 264 : gridder=0;
177 : else{
178 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
179 0 : uvScale, uvOffset,
180 0 : convType);
181 : }
182 264 : convFunc_p.resize();
183 264 : convSupport_p=other.convSupport_p;
184 264 : convSampling_p=other.convSampling_p;
185 264 : convFunc_p=other.convFunc_p;
186 264 : isTiled=other.isTiled;
187 : //lattice=other.lattice;
188 264 : lattice.reset( );
189 264 : tilesize=other.tilesize;
190 264 : arrayLattice.reset( );
191 264 : maxAbsData=other.maxAbsData;
192 264 : centerLoc=other.centerLoc;
193 264 : offsetLoc=other.offsetLoc;
194 264 : padding_p=other.padding_p;
195 264 : usezero_p=other.usezero_p;
196 264 : noPadding_p=other.noPadding_p;
197 264 : machineName_p="GridFT";
198 264 : timemass_p=0.0;
199 264 : timegrid_p=0.0;
200 : // peek = other.peek;
201 : };
202 264 : return *this;
203 : };
204 :
205 : //----------------------------------------------------------------------
206 264 : GridFT::GridFT(const GridFT& other) : FTMachine(), machineName_p("GridFT")
207 : {
208 264 : operator=(other);
209 264 : }
210 :
211 : //-----------------------------------------------------------------------
212 264 : FTMachine* GridFT::cloneFTM(){
213 264 : return new GridFT(*this);
214 : }
215 :
216 : //===============================
217 2141 : Long GridFT::estimateRAM(const CountedPtr<SIImageStore>& imstor){
218 2141 : Long mem=FTMachine::estimateRAM(imstor);
219 2141 : mem=mem*padding_p*padding_p;
220 2141 : return mem;
221 : }
222 :
223 :
224 : //----------------------------------------------------------------------
225 2894 : void GridFT::init() {
226 :
227 2894 : logIO() << LogOrigin("GridFT", "init") << LogIO::NORMAL;
228 :
229 2894 : ok();
230 : // if (peek == NULL)
231 : // {
232 : // // peek = new SynthesisAsyncPeek();
233 : // // peek->reset();
234 : // // peek->startThread();
235 : // }
236 :
237 : /* hardwiring isTiled is false
238 : // Padding is possible only for non-tiled processing
239 : if((padding_p*padding_p*image->shape().product())>cachesize) {
240 : isTiled=true;
241 : nx = image->shape()(0);
242 : ny = image->shape()(1);
243 : npol = image->shape()(2);
244 : nchan = image->shape()(3);
245 : }
246 : else {
247 : */
248 : // We are padding.
249 2894 : isTiled=false;
250 2894 : if(!noPadding_p && padding_p > 1.01){
251 2869 : CompositeNumber cn(uInt(image->shape()(0)*2));
252 2869 : nx = cn.nextLargerEven(Int(padding_p*Float(image->shape()(0))-0.5));
253 2869 : ny = cn.nextLargerEven(Int(padding_p*Float(image->shape()(1))-0.5));
254 2869 : }
255 : else{
256 25 : nx = image->shape()(0);
257 25 : ny = image->shape()(1);
258 : }
259 2894 : npol = image->shape()(2);
260 2894 : nchan = image->shape()(3);
261 : // }
262 :
263 2894 : sumWeight.resize(npol, nchan);
264 :
265 2894 : uvScale.resize(2);
266 2894 : uvScale(0)=(Float(nx)*image->coordinates().increment()(0));
267 2894 : uvScale(1)=(Float(ny)*image->coordinates().increment()(1));
268 2894 : uvOffset.resize(2);
269 2894 : uvOffset(0)=nx/2;
270 2894 : uvOffset(1)=ny/2;
271 :
272 : // Now set up the gridder. The possibilities are BOX and SF
273 2894 : if(gridder) delete gridder; gridder=0;
274 2894 : convType.upcase();
275 5788 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
276 2894 : uvScale, uvOffset,
277 2894 : convType);
278 :
279 2894 : convFunc_p.resize();
280 2894 : convSupport_p=gridder->cSupport()(0);
281 2894 : convSampling_p=gridder->cSampling();
282 2894 : convFunc_p=gridder->cFunction();
283 : // Set up image cache needed for gridding. For BOX-car convolution
284 : // we can use non-overlapped tiles. Otherwise we need to use
285 : // overlapped tiles and additive gridding so that only increments
286 : // to a tile are written.
287 2894 : if(imageCache) delete imageCache; imageCache=0;
288 :
289 2894 : if(isTiled) {
290 0 : Float tileOverlap=0.5;
291 0 : if(convType=="BOX") {
292 0 : tileOverlap=0.0;
293 : }
294 : else {
295 0 : tileOverlap=0.5;
296 0 : tilesize=max(12,tilesize);
297 : }
298 0 : IPosition tileShape=IPosition(4,tilesize,tilesize,npol,nchan);
299 0 : Vector<Float> tileOverlapVec(4);
300 0 : tileOverlapVec=0.0;
301 0 : tileOverlapVec(0)=tileOverlap;
302 0 : tileOverlapVec(1)=tileOverlap;
303 0 : Int tmpCacheVal=static_cast<Int>(cachesize);
304 0 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
305 : tileOverlapVec,
306 0 : (tileOverlap>0.0));
307 :
308 0 : }
309 2894 : }
310 :
311 : // This is nasty, we should use CountedPointers here.
312 5538 : GridFT::~GridFT() {
313 2769 : if(imageCache) delete imageCache; imageCache=0;
314 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
315 2769 : if(gridder) delete gridder; gridder=0;
316 5538 : }
317 :
318 : // Initialize for a transform from the Sky domain. This means that
319 : // we grid-correct, and FFT the image
320 :
321 775 : void GridFT::initializeToVis(ImageInterface<Complex>& iimage,
322 : const vi::VisBuffer2& vb)
323 : {
324 775 : image=&iimage;
325 775 : toVis_p=true;
326 :
327 775 : ok();
328 :
329 775 : init();
330 : // peek->reset();
331 : // Initialize the maps for polarization and channel. These maps
332 : // translate visibility indices into image indices
333 775 : initMaps(vb);
334 : // Need to reset nx, ny for padding
335 : // Padding is possible only for non-tiled processing
336 :
337 :
338 : //cerr << "initialize to vis nx" << nx << " " << ny << " " << npol << " " << nchan << endl;
339 :
340 : //cerr << "image shape " << image->shape() << endl;
341 :
342 : // If we are memory-based then read the image in and create an
343 : // ArrayLattice otherwise just use the PagedImage
344 : /*if(isTiled) {
345 : lattice=std::shared_ptr<Lattice<Complex> >(image, false);
346 : }
347 : else {
348 :
349 : }*/
350 775 : prepGridForDegrid();
351 :
352 : //AlwaysAssert(lattice, AipsError);
353 :
354 :
355 :
356 775 : }
357 :
358 775 : void GridFT::prepGridForDegrid(){
359 775 : IPosition gridShape(4, nx, ny, npol, nchan);
360 775 : griddedData.resize(gridShape);
361 : //griddedData can be a reference of image data...if not using model col
362 : //hence using an undocumented feature of resize that if
363 : //the size is the same as old data it is not changed.
364 : //if(!usePut2_p) griddedData.set(0);
365 775 : griddedData.set(Complex(0.0));
366 :
367 775 : IPosition stride(4, 1);
368 1550 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
369 1550 : IPosition trc(blc+image->shape()-stride);
370 :
371 775 : IPosition start(4, 0);
372 775 : griddedData(blc, trc) = image->getSlice(start, image->shape());
373 :
374 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
375 775 : arrayLattice.reset( new ArrayLattice<Complex>(griddedData) );
376 775 : lattice=arrayLattice;
377 : //logIO() << LogIO::DEBUGGING
378 : // << "Starting grid correction and FFT of image" << LogIO::POST;
379 :
380 : // Do the Grid-correction.
381 : {
382 775 : Vector<Complex> correction(nx);
383 775 : correction=Complex(1.0, 0.0);
384 : // Do the Grid-correction
385 775 : IPosition cursorShape(4, nx, 1, 1, 1);
386 775 : IPosition axisPath(4, 0, 1, 2, 3);
387 775 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
388 775 : LatticeIterator<Complex> lix(*lattice, lsx);
389 547039 : for(lix.reset();!lix.atEnd();lix++) {
390 546264 : gridder->correctX1D(correction, lix.position()(1));
391 546264 : lix.rwVectorCursor()/=correction;
392 : }
393 775 : }
394 775 : image->clearCache();
395 : // Now do the FFT2D in place
396 : //LatticeFFT::cfft2d(*lattice);
397 775 : ft_p.c2cFFT(*lattice);
398 : //logIO() << LogIO::DEBUGGING
399 : // << "Finished grid correction and FFT of image" << LogIO::POST;
400 :
401 775 : }
402 :
403 :
404 773 : void GridFT::finalizeToVis()
405 : {
406 :
407 773 : logIO() << LogOrigin("GridFT", "finalizeToVis") << LogIO::NORMAL;
408 773 : logIO() <<LogIO::NORMAL2<< "Time to degrid " << timedegrid_p <<LogIO::POST;
409 773 : timedegrid_p=0.0;
410 :
411 773 : if(arrayLattice) arrayLattice=nullptr;
412 773 : if(lattice) lattice=nullptr;
413 773 : griddedData.resize();
414 773 : if(isTiled) {
415 :
416 :
417 :
418 0 : AlwaysAssert(imageCache, AipsError);
419 0 : AlwaysAssert(image, AipsError);
420 0 : ostringstream o;
421 0 : imageCache->flush();
422 0 : imageCache->showCacheStatistics(o);
423 0 : logIO() << o.str() << LogIO::POST;
424 0 : }
425 773 : }
426 :
427 :
428 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
429 : // grid.
430 2119 : void GridFT::initializeToSky(ImageInterface<Complex>& iimage,
431 : Matrix<Float>& weight, const vi::VisBuffer2& vb)
432 : {
433 : // image always points to the image
434 2119 : image=&iimage;
435 2119 : toVis_p=false;
436 :
437 2119 : init();
438 :
439 : // Initialize the maps for polarization and channel. These maps
440 : // translate visibility indices into image indices
441 2119 : initMaps(vb);
442 :
443 :
444 :
445 2119 : sumWeight=0.0;
446 2119 : weight.resize(sumWeight.shape());
447 2119 : weight=0.0;
448 :
449 :
450 2119 : IPosition gridShape(4, nx, ny, npol, nchan);
451 2119 : if( !useDoubleGrid_p )
452 : {
453 23 : griddedData.resize(gridShape);
454 23 : griddedData=Complex(0.0);
455 : }
456 : else{
457 :
458 2096 : griddedData2.resize(gridShape);
459 2096 : griddedData2=DComplex(0.0);
460 : }
461 2119 : image->clearCache();
462 : //iimage.get(griddedData, false);
463 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
464 :
465 :
466 : //AlwaysAssert(lattice, AipsError);
467 2119 : }
468 :
469 :
470 :
471 2118 : void GridFT::finalizeToSky()
472 : {
473 : //AlwaysAssert(lattice, AipsError);
474 : // Now we flush the cache and report statistics
475 : // For memory based, we don't write anything out yet.
476 2118 : logIO() << LogOrigin("GridFT", "finalizeToSky") << LogIO::NORMAL;
477 2118 : logIO()<<LogIO::NORMAL2 <<"Time to massage data " << timemass_p << LogIO::POST;
478 2118 : logIO()<< LogIO::NORMAL2 <<"Time to grid data " << timegrid_p << LogIO::POST;
479 2118 : timemass_p=0.0;
480 2118 : timegrid_p=0.0;
481 2118 : if(isTiled) {
482 :
483 :
484 0 : AlwaysAssert(image, AipsError);
485 0 : AlwaysAssert(imageCache, AipsError);
486 0 : imageCache->flush();
487 0 : ostringstream o;
488 0 : imageCache->showCacheStatistics(o);
489 0 : logIO() << o.str() << LogIO::POST;
490 0 : }
491 : // peek->terminate();
492 : // peek->reset();
493 2118 : }
494 :
495 :
496 :
497 0 : Array<Complex>* GridFT::getDataPointer(const IPosition& centerLoc2D,
498 : Bool readonly) {
499 : Array<Complex>* result;
500 : // Is tiled: get tiles and set up offsets
501 0 : centerLoc(0)=centerLoc2D(0);
502 0 : centerLoc(1)=centerLoc2D(1);
503 0 : result=&imageCache->tile(offsetLoc,centerLoc, readonly);
504 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
505 0 : return result;
506 : }
507 :
508 : #define NEED_UNDERSCORES
509 : #if defined(NEED_UNDERSCORES)
510 : #define ggrid ggrid_
511 : #define dgrid dgrid_
512 : #define ggrids ggrids_
513 : #define sectggridd sectggridd_
514 : #define sectggrids sectggrids_
515 : #define sectdgrid sectdgrid_
516 : #define locuvw locuvw_
517 : #endif
518 :
519 : extern "C" {
520 : void locuvw(const Double*, const Double*, const Double*, const Int*, const Double*, const Double*, const Int*,
521 : Int*, Int*, Complex*, const Int*, const Int*, const Double*);
522 : void ggrid(Double*,
523 : Double*,
524 : const Complex*,
525 : Int*,
526 : Int*,
527 : Int*,
528 : const Int*,
529 : const Int*,
530 : const Float*,
531 : Int*,
532 : Int*,
533 : Double*,
534 : Double*,
535 : DComplex*,
536 : Int*,
537 : Int*,
538 : Int *,
539 : Int *,
540 : const Double*,
541 : const Double*,
542 : Int*,
543 : Int*,
544 : Double*,
545 : Int*,
546 : Int*,
547 : Double*);
548 : void sectggridd(const Complex*,
549 : const Int*,
550 : const Int*,
551 : const Int*,
552 : const Int*,
553 : const Int*,
554 : const Float*,
555 : const Int*,
556 : DComplex*,
557 : const Int*,
558 : const Int*,
559 : const Int *,
560 : const Int *,
561 : //support
562 : const Int*,
563 : const Int*,
564 : const Double*,
565 : const Int*,
566 : const Int*,
567 : Double*,
568 : const Int*,
569 : const Int*,
570 : const Int*,
571 : const Int*,
572 : const Int*,
573 : const Int*,
574 : //loc, off, phasor
575 : const Int*,
576 : const Int*,
577 : const Complex*);
578 : //single precision gridder
579 : void sectggrids(const Complex*,
580 : const Int*,
581 : const Int*,
582 : const Int*,
583 : const Int*,
584 : const Int*,
585 : const Float*,
586 : const Int*,
587 : Complex*,
588 : const Int*,
589 : const Int*,
590 : const Int *,
591 : const Int *,
592 : const Int*,
593 : const Int*,
594 : const Double*,
595 : const Int*,
596 : const Int*,
597 : Double*,
598 : const Int*,
599 : const Int*,
600 : const Int*,
601 : const Int*,
602 : const Int*,
603 : const Int*,
604 : //loc, off, phasor
605 : const Int*,
606 : const Int*,
607 : const Complex*);
608 : void ggrids(Double*,
609 : Double*,
610 : const Complex*,
611 : Int*,
612 : Int*,
613 : Int*,
614 : const Int*,
615 : const Int*,
616 : const Float*,
617 : Int*,
618 : Int*,
619 : Double*,
620 : Double*,
621 : Complex*,
622 : Int*,
623 : Int*,
624 : Int *,
625 : Int *,
626 : const Double*,
627 : const Double*,
628 : Int*,
629 : Int*,
630 : Double*,
631 : Int*,
632 : Int*,
633 : Double*);
634 :
635 : void dgrid(Double*,
636 : Double*,
637 : Complex*,
638 : Int*,
639 : Int*,
640 : const Int*,
641 : const Int*,
642 : Int*,
643 : Int*,
644 : Double*,
645 : Double*,
646 : const Complex*,
647 : Int*,
648 : Int*,
649 : Int *,
650 : Int *,
651 : const Double*,
652 : const Double*,
653 : Int*,
654 : Int*,
655 : Double*,
656 : Int*,
657 : Int*);
658 :
659 : void sectdgrid(Complex*,
660 : const Int*,
661 : const Int*,
662 : const Int*,
663 : const Int*,
664 : const Int*,
665 : const Complex*,
666 : const Int*,
667 : const Int*,
668 : const Int *,
669 : const Int *,
670 : //support
671 : const Int*,
672 : const Int*,
673 : const Double*,
674 : const Int*,
675 : const Int*,
676 : //rbeg, rend, loc, off, phasor
677 : const Int*,
678 : const Int*,
679 : const Int*,
680 : const Int*,
681 : const Complex*);
682 : }
683 509336 : void GridFT::put(const vi::VisBuffer2& vb, Int row, Bool dopsf,
684 : FTMachine::Type type)
685 : {
686 :
687 509336 : gridOk(convSupport_p);
688 : // peek->setVBPtr(&vb);
689 :
690 : //Check if ms has changed then cache new spw and chan selection
691 : //if(vb.isNewMs())
692 : // matchAllSpwChans(vb);
693 :
694 : //Here we redo the match or use previous match
695 :
696 : //Channel matching for the actual spectral window of buffer
697 : //if(doConversion_p[vb.spectralWindows()[0]]){
698 509335 : matchChannel(vb);
699 : //}
700 : //else{
701 : // chanMap.resize();
702 : // chanMap=multiChanMap_p[vb.spectralWindows()[0]];
703 : //}
704 :
705 : //No point in reading data if its not matching in frequency
706 509335 : if(max(chanMap)==-1)
707 : {
708 : // peek->reset();
709 136 : return;
710 : }
711 :
712 509199 : Timer tim;
713 509199 : tim.mark();
714 :
715 : //const Matrix<Float> *imagingweight;
716 : //imagingweight=&(vb.imagingWeight());
717 509199 : Matrix<Float> imagingweight;
718 509199 : getImagingWeight(imagingweight, vb);
719 :
720 509199 : if(dopsf) {type=FTMachine::PSF;}
721 :
722 509199 : Cube<Complex> data;
723 : //Fortran gridder need the flag as ints
724 509199 : Cube<Int> flags;
725 509199 : Matrix<Float> elWeight;
726 509199 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
727 :
728 :
729 : Bool iswgtCopy;
730 : const Float *wgtStorage;
731 509199 : wgtStorage=elWeight.getStorage(iswgtCopy);
732 :
733 : Bool isCopy;
734 : const Complex *datStorage;
735 509199 : if(!dopsf)
736 311538 : datStorage=data.getStorage(isCopy);
737 : else
738 197661 : datStorage=0;
739 : // If row is -1 then we pass through all rows
740 : Int startRow, endRow, nRow;
741 509199 : if (row==-1) {
742 509199 : nRow=vb.nRows();
743 509199 : startRow=0;
744 509199 : endRow=nRow-1;
745 : } else {
746 0 : nRow=1;
747 0 : startRow=row;
748 0 : endRow=row;
749 : }
750 :
751 : //cerr << "nRow " << nRow << endl;
752 :
753 :
754 : // Get the uvws in a form that Fortran can use and do that
755 : // necessary phase rotation. On a Pentium Pro 200 MHz
756 : // when null, this step takes about 50us per uvw point. This
757 : // is just barely noticeable for Stokes I continuum and
758 : // irrelevant for other cases.
759 :
760 : Bool del;
761 : // IPosition s(flags.shape());
762 509199 : const IPosition &fs=flags.shape();
763 509199 : std::vector<Int> s(fs.begin(),fs.end());
764 509199 : Int nvispol=s[0];
765 509199 : Int nvischan=s[1];
766 509199 : Int nvisrow=s[2];
767 509199 : Matrix<Double> uvw(negateUV(vb));
768 :
769 509199 : Vector<Double> dphase(vb.nRows());
770 509199 : Cube<Int> loc(2, nvischan, nRow);
771 509199 : Cube<Int> off(2, nvischan, nRow);
772 509199 : Matrix<Complex> phasor(nvischan, nRow);
773 509199 : Int csamp=convSampling_p;
774 509199 : dphase=0.0;
775 :
776 509199 : timemass_p +=tim.real();
777 509199 : tim.mark();
778 509199 : rotateUVW(uvw, dphase, vb);
779 509199 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
780 : Bool delphase;
781 509199 : Complex * phasorstor=phasor.getStorage(delphase);
782 509199 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
783 509199 : const Double * scalestor=uvScale.getStorage(del);
784 509199 : const Double * offsetstor=uvOffset.getStorage(del);
785 509199 : const Double* uvstor= uvw.getStorage(del);
786 509199 : Int * locstor=loc.getStorage(del);
787 509199 : Int * offstor=off.getStorage(del);
788 509199 : const Double *dpstor=dphase.getStorage(del);
789 : //Vector<Double> f1=interpVisFreq_p.copy();
790 509199 : Int nvchan=nvischan;
791 : Int irow;
792 509199 : Double cinv=Double(1.0)/C::c;
793 509199 : Int dow=0;
794 509199 : Int nth=1;
795 : #ifdef _OPENMP
796 509199 : if(numthreads_p >0){
797 0 : nth=min(numthreads_p, omp_get_max_threads());
798 : }
799 : else{
800 509199 : nth= omp_get_max_threads();
801 : }
802 : //nth=min(4,nth);
803 : #endif
804 :
805 :
806 509199 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvchan, scalestor, offsetstor, csamp, phasorstor, uvstor, locstor, offstor, dpstor, cinv, dow) shared(startRow, endRow) num_threads(nth)
807 : {
808 : #pragma omp for schedule(dynamic)
809 : for (irow=startRow; irow<=endRow; ++irow){
810 : //locateuvw(uvstor,dpstor, visfreqstor, nvchan, scalestor, offsetstor, csamp,
811 : // locstor,
812 : // offstor, phasorstor, irow);
813 : locuvw(uvstor, dpstor, visfreqstor, &nvchan, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
814 : }
815 :
816 : }//end pragma parallel
817 : // Take care of translation of Bools to Integer
818 509199 : Int idopsf=0;
819 509199 : if(dopsf) idopsf=1;
820 :
821 : //////TESTOO
822 : //ofstream myfile;
823 : //myfile.open ("putLoc.txt", ios::out | ios::app | ios::ate );
824 : //myfile << vb.rowIds()(0) << " uv " << uvw.column(0) << " loc " << loc(0,0,0) << ", " << loc(1,0,0) << "\n" << endl;
825 : //myfile.close();
826 : ///////////////
827 :
828 :
829 :
830 509199 : Vector<Int> rowFlags(vb.nRows());
831 509199 : rowFlags=0;
832 509199 : rowFlags(vb.flagRow())=true;
833 509199 : if(!usezero_p) {
834 172382192 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
835 171878527 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
836 : }
837 : }
838 :
839 : /////////////Some extra stuff for openmp
840 :
841 : Int ixsub, iysub, icounter;
842 509199 : Int csupp=convSupport_p;
843 :
844 509199 : const Double * convfuncstor=(convFunc_p).getStorage(del);
845 :
846 : // cerr <<"Poffset " << min(off) << " " << max(off) << " length " << gridder->cFunction().shape() << endl;
847 :
848 :
849 :
850 509199 : ixsub=1;
851 509199 : iysub=1;
852 509199 : if (nth >4){
853 244278 : ixsub=8;
854 244278 : iysub=8;
855 : }
856 264921 : else if(nth >1) {
857 264921 : ixsub=2;
858 264921 : iysub=2;
859 : }
860 :
861 :
862 509199 : Int rbeg=startRow+1;
863 509199 : Int rend=endRow+1;
864 509199 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
865 17202675 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
866 16693476 : sumwgt[icounter].resize(sumWeight.shape());
867 16693476 : sumwgt[icounter].set(0.0);
868 : }
869 509199 : if(doneThreadPartition_p < 0){
870 1175 : xsect_p.resize(ixsub*iysub);
871 1175 : ysect_p.resize(ixsub*iysub);
872 1175 : nxsect_p.resize(ixsub*iysub);
873 1175 : nysect_p.resize(ixsub*iysub);
874 36475 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
875 35300 : findGridSector(nx, ny, ixsub, iysub, 0, 0, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
876 : }
877 : }
878 509199 : Vector<Int> xsect, ysect, nxsect, nysect;
879 509199 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
880 509199 : const Int* pmapstor=polMap.getStorage(del);
881 509199 : const Int *cmapstor=chanMap.getStorage(del);
882 509199 : Int nc=nchan;
883 509199 : Int np=npol;
884 509199 : Int nxp=nx;
885 509199 : Int nyp=ny;
886 509199 : const Int * flagstor=flags.getStorage(del);
887 509199 : const Int * rowflagstor=rowFlags.getStorage(del);
888 : ////////////////////////
889 :
890 : Bool gridcopy;
891 509199 : if(useDoubleGrid_p){
892 503265 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
893 503265 : #pragma omp parallel default(none) private(icounter, del) firstprivate(idopsf, datStorage, wgtStorage, flagstor, rowflagstor, convfuncstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc,ixsub, iysub, rend, rbeg, csamp, csupp, nvispol, nvischan, nvisrow, phasorstor, locstor, offstor, xsect, ysect, nxsect, nysect) shared(sumwgt) num_threads(nth)
894 :
895 : {
896 : //cerr << "numthreads " << omp_get_num_threads() << endl;
897 : #pragma omp for schedule(dynamic)
898 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
899 : //cerr << "thread id " << omp_get_thread_num() << endl;
900 : Int x0=xsect(icounter);
901 : Int y0=ysect(icounter);
902 : Int nxsub=nxsect(icounter);
903 : Int nysub=nysect(icounter);
904 :
905 : sectggridd(datStorage,
906 : &nvispol,
907 : &nvischan,
908 : &idopsf,
909 : flagstor,
910 : rowflagstor,
911 : wgtStorage,
912 : &nvisrow,
913 : gridstor,
914 : &nxp,
915 : &nyp,
916 : &np,
917 : &nc,
918 : &csupp,
919 : &csamp,
920 : convfuncstor,
921 : cmapstor,
922 : pmapstor,
923 : (sumwgt[icounter]).getStorage(del),
924 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
925 : phasorstor);
926 : }//end for
927 : }// end pragma parallel
928 17120205 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
929 16616940 : sumWeight=sumWeight+sumwgt[icounter];
930 : }
931 : //phasor.putStorage(phasorstor, delphase);
932 503265 : griddedData2.putStorage(gridstor, gridcopy);
933 503265 : if(dopsf && (nth >4))
934 99198 : tweakGridSector(nx, ny, ixsub, iysub);
935 : }
936 : else{
937 5934 : Complex *gridstor=griddedData.getStorage(gridcopy);
938 5934 : #pragma omp parallel default(none) private(icounter, del) firstprivate(idopsf, datStorage, wgtStorage, flagstor, rowflagstor, convfuncstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc,ixsub, iysub, rend, rbeg, csamp, csupp, nvispol, nvischan, nvisrow, phasorstor, locstor, offstor, xsect, ysect, nxsect, nysect) shared(sumwgt) num_threads(ixsub*iysub)
939 : {
940 : //cerr << "numthreads " << omp_get_num_threads() << endl;
941 : #pragma omp for schedule(dynamic)
942 :
943 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
944 : //cerr << "thread id " << omp_get_thread_num() << endl;
945 : Int x0=xsect(icounter);
946 : Int y0=ysect(icounter);
947 : Int nxsub=nxsect(icounter);
948 : Int nysub=nysect(icounter);
949 :
950 :
951 :
952 : //cerr << "x0 " << x0 << " y0 " << y0 << " nxsub " << nxsub << " nysub " << nysub << endl;
953 : sectggrids(datStorage,
954 : &nvispol,
955 : &nvischan,
956 : &idopsf,
957 : flagstor,
958 : rowflagstor,
959 : wgtStorage,
960 : &nvisrow,
961 : gridstor,
962 : &nxp,
963 : &nyp,
964 : &np,
965 : &nc,
966 : &csupp,
967 : &csamp,
968 : convfuncstor,
969 : cmapstor,
970 : pmapstor,
971 : (sumwgt[icounter]).getStorage(del),
972 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
973 : phasorstor);
974 : }//end for
975 : }// end pragma parallel
976 82470 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
977 76536 : sumWeight=sumWeight+sumwgt[icounter];
978 : }
979 :
980 5934 : griddedData.putStorage(gridstor, gridcopy);
981 5934 : if(dopsf && (nth > 4))
982 480 : tweakGridSector(nx, ny, ixsub, iysub);
983 : }
984 : // cerr << "sunweight " << sumWeight << endl;
985 :
986 509199 : timegrid_p+=tim.real();
987 :
988 509199 : if(!dopsf)
989 311538 : data.freeStorage(datStorage, isCopy);
990 509199 : elWeight.freeStorage(wgtStorage,iswgtCopy);
991 :
992 : // peek->reset();
993 509199 : }
994 :
995 21 : void GridFT::modifyConvFunc(const Vector<Double>& convFunc, Int convSupport, Int convSampling){
996 21 : convFunc_p.resize();
997 21 : convFunc_p=convFunc;
998 21 : convSupport_p=convSupport;
999 21 : convSampling_p=convSampling;
1000 :
1001 21 : }
1002 :
1003 159984 : void GridFT::get(vi::VisBuffer2& vb, Int row)
1004 : {
1005 :
1006 159984 : gridOk(convSupport_p);
1007 : // If row is -1 then we pass through all rows
1008 : Int startRow, endRow, nRow;
1009 159984 : if (row < 0) {
1010 159984 : nRow=vb.nRows();
1011 159984 : startRow=0;
1012 159984 : endRow=nRow-1;
1013 : //unnecessary zeroing
1014 : //vb.modelVisCube()=Complex(0.0,0.0);
1015 : } else {
1016 0 : nRow=1;
1017 0 : startRow=row;
1018 0 : endRow=row;
1019 : //unnecessary zeroing
1020 : //vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
1021 : }
1022 :
1023 :
1024 : ///Channel matching for the actual spectral window of buffer
1025 159984 : matchChannel(vb);
1026 :
1027 :
1028 : //cerr << "GETchanMap " << chanMap << endl;
1029 : //No point in reading data if its not matching in frequency
1030 159984 : if(max(chanMap)==-1)
1031 58 : return;
1032 :
1033 :
1034 :
1035 : // Get the uvws in a form that Fortran can use
1036 159926 : Matrix<Double> uvw(negateUV(vb));
1037 159926 : Vector<Double> dphase(vb.nRows());
1038 159926 : dphase=0.0;
1039 159926 : rotateUVW(uvw, dphase, vb);
1040 159926 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1041 :
1042 :
1043 :
1044 : //Check if ms has changed then cache new spw and chan selection
1045 : //if(vb.isNewMs())
1046 : // matchAllSpwChans(vb);
1047 :
1048 :
1049 : //Here we redo the match or use previous match
1050 :
1051 :
1052 159926 : Cube<Complex> data;
1053 159926 : Cube<Int> flags;
1054 159926 : getInterpolateArrays(vb, data, flags);
1055 :
1056 : //cerr << "max min data " << max(griddedData) << " " << min(griddedData) << " flags " << max(flags) << " " << min(flags) << endl;
1057 :
1058 : // IPosition s(data.shape());
1059 159926 : Int nvp=data.shape()(0);
1060 159926 : Int nvc=data.shape()(1);
1061 159926 : Int nvisrow=data.shape()(2);
1062 :
1063 :
1064 : //cerr << "get flags " << min(flags) << " " << max(flags) << endl;
1065 : Complex *datStorage;
1066 : Bool isCopy;
1067 159926 : datStorage=data.getStorage(isCopy);
1068 :
1069 : ///
1070 159926 : Cube<Int> loc(2, nvc, nvisrow);
1071 159926 : Cube<Int> off(2, nvc, nRow);
1072 159926 : Matrix<Complex> phasor(nvc, nRow);
1073 159926 : Int csamp=convSampling_p;
1074 : Bool delphase;
1075 : Bool del;
1076 159926 : Complex * phasorstor=phasor.getStorage(delphase);
1077 159926 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1078 159926 : const Double * scalestor=uvScale.getStorage(del);
1079 159926 : const Double * offsetstor=uvOffset.getStorage(del);
1080 159926 : const Double* uvstor= uvw.getStorage(del);
1081 159926 : Int * locstor=loc.getStorage(del);
1082 159926 : Int * offstor=off.getStorage(del);
1083 159926 : const Double *dpstor=dphase.getStorage(del);
1084 : //Vector<Double> f1=interpVisFreq_p.copy();
1085 159926 : Int nvchan=nvc;
1086 : Int irow;
1087 159926 : Double cinv=Double(1.0)/C::c;
1088 159926 : Int dow=0;
1089 159926 : Int nth=1;
1090 : #ifdef _OPENMP
1091 159926 : if(numthreads_p >0){
1092 0 : nth=min(numthreads_p, omp_get_max_threads());
1093 : }
1094 : else{
1095 159926 : nth= omp_get_max_threads();
1096 : }
1097 : //nth=min(4,nth);
1098 : #endif
1099 :
1100 :
1101 159926 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvchan, scalestor, offsetstor, csamp, phasorstor, uvstor, locstor, offstor, dpstor, cinv, dow) shared(startRow, endRow) num_threads(nth)
1102 :
1103 : {
1104 : #pragma omp for schedule(dynamic)
1105 : for (irow=startRow; irow<=endRow; ++irow){
1106 : //locateuvw(uvstor,dpstor, visfreqstor, nvchan, scalestor, offsetstor, csamp,
1107 : // locstor,
1108 : // offstor, phasorstor, irow);
1109 :
1110 : locuvw(uvstor, dpstor, visfreqstor, &nvchan, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
1111 : }
1112 :
1113 : }//end pragma parallel
1114 :
1115 159926 : Int rbeg=startRow+1;
1116 159926 : Int rend=endRow+1;
1117 :
1118 :
1119 159926 : Vector<Int> rowFlags(vb.nRows());
1120 159926 : rowFlags=0;
1121 159926 : rowFlags(vb.flagRow())=true;
1122 : //cerr << "rowFlags " << rowFlags << endl;
1123 159926 : if(!usezero_p) {
1124 55143110 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1125 54983184 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
1126 : }
1127 : }
1128 :
1129 :
1130 : //cerr <<"offset " << min(off) << " " <<max(off) << " length " << gridder->cFunction().shape() << endl;
1131 :
1132 : {
1133 : Bool delgrid;
1134 159926 : const Complex* gridstor=griddedData.getStorage(delgrid);
1135 159926 : const Double * convfuncstor=(convFunc_p).getStorage(del);
1136 :
1137 159926 : const Int* pmapstor=polMap.getStorage(del);
1138 159926 : const Int *cmapstor=chanMap.getStorage(del);
1139 159926 : Int nc=nchan;
1140 159926 : Int np=npol;
1141 159926 : Int nxp=nx;
1142 159926 : Int nyp=ny;
1143 159926 : Int csupp=convSupport_p;
1144 159926 : const Int * flagstor=flags.getStorage(del);
1145 159926 : const Int * rowflagstor=rowFlags.getStorage(del);
1146 :
1147 :
1148 159926 : Int npart=nth;
1149 :
1150 :
1151 159926 : Int ix=0;
1152 159926 : #pragma omp parallel default(none) private(ix, rbeg, rend) firstprivate(datStorage, flagstor, rowflagstor, convfuncstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc, csamp, csupp, nvp, nvc, nvisrow, phasorstor, locstor, offstor) shared(npart) num_threads(npart)
1153 : {
1154 : #pragma omp for schedule(dynamic)
1155 : for (ix=0; ix< npart; ++ix){
1156 : rbeg=ix*(nvisrow/npart)+1;
1157 : rend=(ix != (npart-1)) ? (rbeg+(nvisrow/npart)-1) : (rbeg+(nvisrow/npart)+nvisrow%npart-1) ;
1158 : //cerr << "rbeg " << rbeg << " rend " << rend << " " << nvisrow << endl;
1159 :
1160 : sectdgrid(datStorage,
1161 : &nvp,
1162 : &nvc,
1163 : flagstor,
1164 : rowflagstor,
1165 : &nvisrow,
1166 : gridstor,
1167 : &nxp,
1168 : &nyp,
1169 : &np,
1170 : &nc,
1171 : &csupp,
1172 : &csamp,
1173 : convfuncstor,
1174 : cmapstor,
1175 : pmapstor,
1176 : &rbeg, &rend, locstor, offstor, phasorstor);
1177 : }//end pragma parallel
1178 : }
1179 159926 : data.putStorage(datStorage, isCopy);
1180 159926 : griddedData.freeStorage(gridstor, delgrid);
1181 : //cerr << "Get min max " << min(data) << " " << max(data) << endl;
1182 : }
1183 159926 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1184 :
1185 159926 : }
1186 :
1187 :
1188 :
1189 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1190 : // return the resulting image
1191 2097 : ImageInterface<Complex>& GridFT::getImage(Matrix<Float>& weights, Bool normalize)
1192 : {
1193 : //AlwaysAssert(lattice, AipsError);
1194 2097 : AlwaysAssert(gridder, AipsError);
1195 2097 : AlwaysAssert(image, AipsError);
1196 2097 : logIO() << LogOrigin("GridFT", "getImage") << LogIO::NORMAL;
1197 :
1198 2097 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
1199 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
1200 : }
1201 2097 : weights.resize(sumWeight.shape());
1202 :
1203 2097 : convertArray(weights, sumWeight);
1204 : // If the weights are all zero then we cannot normalize
1205 : // otherwise we don't care.
1206 2097 : if(normalize&&max(weights)==0.0) {
1207 : logIO() << LogIO::SEVERE << "No useful data in GridFT: weights all zero"
1208 0 : << LogIO::POST;
1209 : }
1210 : else {
1211 :
1212 : //const IPosition latticeShape = lattice->shape();
1213 :
1214 : //logIO() << LogIO::DEBUGGING
1215 : // << "Starting FFT and scaling of image" << LogIO::POST;
1216 :
1217 :
1218 :
1219 : // if(useDoubleGrid_p){
1220 : // convertArray(griddedData, griddedData2);
1221 : // //Don't need the double-prec grid anymore...
1222 : // griddedData2.resize();
1223 : // }
1224 :
1225 : // x and y transforms
1226 : // LatticeFFT::cfft2d(*lattice,false);
1227 : //
1228 : // Retain the double precision grid for FFT as well. Convert it
1229 : // to single precision just after (since images are still single
1230 : // precision).
1231 : //
1232 2097 : if(useDoubleGrid_p)
1233 : {
1234 2095 : ArrayLattice<DComplex> darrayLattice(griddedData2);
1235 : //LatticeFFT::cfft2d(darrayLattice,false);
1236 2095 : ft_p.c2cFFT(darrayLattice, False);
1237 2095 : griddedData.resize(griddedData2.shape());
1238 2095 : convertArray(griddedData, griddedData2);
1239 :
1240 2095 : SynthesisUtilMethods::getResource("mem peak in getImage");
1241 :
1242 : //Don't need the double-prec grid anymore...
1243 2095 : griddedData2.resize();
1244 2095 : arrayLattice.reset( new ArrayLattice<Complex>(griddedData) );
1245 2095 : lattice=arrayLattice;
1246 2095 : }
1247 : else{
1248 2 : arrayLattice.reset( new ArrayLattice<Complex>(griddedData) );
1249 2 : lattice=arrayLattice;
1250 :
1251 : //LatticeFFT::cfft2d(*lattice,false);
1252 2 : ft_p.c2cFFT(*lattice, False);
1253 : }
1254 :
1255 :
1256 :
1257 : {
1258 2097 : Int inx = lattice->shape()(0);
1259 2097 : Int iny = lattice->shape()(1);
1260 2097 : Vector<Complex> correction(inx);
1261 2097 : correction=Complex(1.0, 0.0);
1262 : // Do the Grid-correction
1263 2097 : IPosition cursorShape(4, inx, 1, 1, 1);
1264 2097 : IPosition axisPath(4, 0, 1, 2, 3);
1265 2097 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1266 2097 : LatticeIterator<Complex> lix(*lattice, lsx);
1267 1547545 : for(lix.reset();!lix.atEnd();lix++) {
1268 1545448 : Int pol=lix.position()(2);
1269 1545448 : Int chan=lix.position()(3);
1270 1545448 : if(weights(pol, chan)!=0.0) {
1271 1499440 : gridder->correctX1D(correction, lix.position()(1));
1272 1499440 : lix.rwVectorCursor()/=correction;
1273 1499440 : if(normalize) {
1274 200 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
1275 200 : lix.rwCursor()*=rnorm;
1276 : }
1277 : else {
1278 1499240 : Complex rnorm(Float(inx)*Float(iny));
1279 1499240 : lix.rwCursor()*=rnorm;
1280 : }
1281 : }
1282 : else {
1283 46008 : lix.woCursor()=0.0;
1284 : }
1285 : }
1286 2097 : }
1287 2097 : LatticeLocker lock1 (*(image), FileLocker::Write);
1288 2097 : if(!isTiled) {
1289 : // Check the section from the image BEFORE converting to a lattice
1290 4194 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1291 2097 : IPosition stride(4, 1);
1292 4194 : IPosition trc(blc+image->shape()-stride);
1293 : // Do the copy
1294 2097 : IPosition start(4, 0);
1295 :
1296 2097 : image->put(griddedData(blc, trc));
1297 2097 : }
1298 2097 : }
1299 2097 : image->flush();
1300 2097 : image->clearCache();
1301 :
1302 2097 : if(arrayLattice) arrayLattice=nullptr;
1303 2097 : if(lattice) lattice=nullptr;
1304 2097 : griddedData.resize();
1305 2097 : return *image;
1306 : }
1307 :
1308 : // Get weight image
1309 0 : void GridFT::getWeightImage(ImageInterface<Float>& weightImage, Matrix<Float>& weights)
1310 : {
1311 :
1312 0 : logIO() << LogOrigin("GridFT", "getWeightImage") << LogIO::NORMAL;
1313 :
1314 0 : weights.resize(sumWeight.shape());
1315 0 : convertArray(weights,sumWeight);
1316 :
1317 0 : const IPosition latticeShape = weightImage.shape();
1318 :
1319 0 : Int nx=latticeShape(0);
1320 0 : Int ny=latticeShape(1);
1321 :
1322 0 : IPosition loc(2, 0);
1323 0 : IPosition cursorShape(4, nx, ny, 1, 1);
1324 0 : IPosition axisPath(4, 0, 1, 2, 3);
1325 0 : LatticeStepper lsx(latticeShape, cursorShape, axisPath);
1326 0 : LatticeIterator<Float> lix(weightImage, lsx);
1327 0 : for(lix.reset();!lix.atEnd();lix++) {
1328 0 : Int pol=lix.position()(2);
1329 0 : Int chan=lix.position()(3);
1330 0 : lix.rwCursor()=weights(pol,chan);
1331 : }
1332 0 : }
1333 :
1334 42 : Bool GridFT::toRecord(String& error,
1335 : RecordInterface& outRec, Bool withImage, const String diskimage)
1336 : {
1337 :
1338 :
1339 :
1340 : // Save the current GridFT object to an output state record
1341 42 : Bool retval = true;
1342 42 : Float elpadd=padding_p;
1343 42 : if(toVis_p && withImage)
1344 19 : elpadd=1.0;
1345 : //save the base class variables
1346 42 : if(!FTMachine::toRecord(error, outRec, withImage, diskimage))
1347 0 : return false;
1348 :
1349 : //a call to init will redo imagecache and gridder
1350 : // so no need to save these
1351 :
1352 42 : outRec.define("cachesize", Int64(cachesize));
1353 42 : outRec.define("tilesize", tilesize);
1354 42 : outRec.define("convtype", convType);
1355 42 : outRec.define("uvscale", uvScale);
1356 42 : outRec.define("uvoffset", uvOffset);
1357 42 : outRec.define("istiled", isTiled);
1358 42 : outRec.define("maxabsdata", maxAbsData);
1359 42 : Vector<Int> center_loc(4), offset_loc(4);
1360 210 : for (Int k=0; k<4 ; k++){
1361 168 : center_loc(k)=centerLoc(k);
1362 168 : offset_loc(k)=offsetLoc(k);
1363 : }
1364 42 : outRec.define("centerloc", center_loc);
1365 42 : outRec.define("offsetloc", offset_loc);
1366 42 : outRec.define("padding", elpadd);
1367 42 : outRec.define("usezero", usezero_p);
1368 42 : outRec.define("nopadding", noPadding_p);
1369 42 : return retval;
1370 42 : }
1371 :
1372 0 : Bool GridFT::fromRecord(String& error,
1373 : const RecordInterface& inRec)
1374 : {
1375 0 : Bool retval = true;
1376 0 : if(!FTMachine::fromRecord(error, inRec))
1377 0 : return false;
1378 0 : gridder=0; imageCache=0; lattice.reset( ); arrayLattice.reset( );
1379 : //For some reason int64 seems to be getting lost...
1380 : //this is not an important parameter...so filing a jira
1381 0 : if(inRec.isDefined("cachesize"))
1382 0 : cachesize=inRec.asInt64("cachesize");
1383 : else
1384 0 : cachesize=1000000;
1385 0 : inRec.get("tilesize", tilesize);
1386 0 : inRec.get("convtype", convType);
1387 0 : inRec.get("uvscale", uvScale);
1388 0 : inRec.get("uvoffset", uvOffset);
1389 0 : inRec.get("istiled", isTiled);
1390 0 : inRec.get("maxabsdata", maxAbsData);
1391 0 : Vector<Int> center_loc(4), offset_loc(4);
1392 0 : inRec.get("centerloc", center_loc);
1393 0 : inRec.get("offsetloc", offset_loc);
1394 0 : uInt ndim4 = 4;
1395 0 : centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2),
1396 0 : center_loc(3));
1397 0 : offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2),
1398 0 : offset_loc(3));
1399 0 : inRec.get("padding", padding_p);
1400 0 : inRec.get("usezero", usezero_p);
1401 0 : inRec.get("nopadding", noPadding_p);
1402 :
1403 0 : machineName_p="GridFT";
1404 : ///setup some of the parameters if there is an image
1405 0 : if(inRec.isDefined("image") || inRec.isDefined("diskimage"))
1406 0 : init();
1407 : /*if(!cmplxImage_p.null()){
1408 : //FTMachine::fromRecord would have recovered the image
1409 : // Might be changing the shape of sumWeight
1410 :
1411 : if(isTiled) {
1412 : lattice=std::shared_ptr<Lattice<Complex> >(image, false);
1413 : }
1414 : else {
1415 : // Make the grid the correct shape and turn it into an array lattice
1416 : // Check the section from the image BEFORE converting to a lattice
1417 : if(!toVis_p){
1418 : IPosition gridShape(4, nx, ny, npol, nchan);
1419 : griddedData.resize(gridShape);
1420 : griddedData=Complex(0.0);
1421 : }
1422 : else{
1423 : prepGridForDegrid();
1424 : }
1425 : }
1426 : };*/
1427 0 : return retval;
1428 0 : }
1429 :
1430 1342319 : void GridFT::ok() {
1431 1342319 : AlwaysAssert(image, AipsError);
1432 1342319 : }
1433 :
1434 : // Make a plain straightforward honest-to-God image. This returns
1435 : // a complex image, without conversion to Stokes. The representation
1436 : // is that required for the visibilities.
1437 : //----------------------------------------------------------------------
1438 0 : void GridFT::makeImage(FTMachine::Type type,
1439 : vi::VisibilityIterator2& vi,
1440 : ImageInterface<Complex>& theImage,
1441 : Matrix<Float>& weight) {
1442 :
1443 :
1444 0 : logIO() << LogOrigin("GridFT", "makeImage") << LogIO::NORMAL;
1445 :
1446 0 : if(type==FTMachine::COVERAGE) {
1447 0 : logIO() << "Type COVERAGE not defined for Fourier transforms" << LogIO::EXCEPTION;
1448 : }
1449 :
1450 :
1451 :
1452 :
1453 : // Loop over all visibilities and pixels
1454 0 : vi::VisBuffer2* vb=vi.getVisBuffer();
1455 :
1456 : // Initialize put (i.e. transform to Sky) for this model
1457 0 : vi.origin();
1458 :
1459 0 : if(vb->polarizationFrame()==MSIter::Linear) {
1460 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
1461 : }
1462 : else {
1463 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
1464 : }
1465 :
1466 0 : initializeToSky(theImage,weight,*vb);
1467 :
1468 : // Loop over the visibilities, putting VisBuffers
1469 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1470 0 : for (vi.origin(); vi.more(); vi.next()) {
1471 :
1472 0 : switch(type) {
1473 0 : case FTMachine::RESIDUAL:
1474 0 : vb->setVisCube(vb->visCubeCorrected()-vb->visCubeModel());
1475 0 : put(*vb, -1, false);
1476 0 : break;
1477 0 : case FTMachine::MODEL:
1478 0 : vb->setVisCube(vb->visCubeModel());
1479 0 : put(*vb, -1, false);
1480 0 : break;
1481 0 : case FTMachine::CORRECTED:
1482 0 : vb->setVisCube(vb->visCubeCorrected());
1483 0 : put(*vb, -1, false);
1484 0 : break;
1485 0 : case FTMachine::PSF:
1486 0 : vb->setVisCube(Complex(1.0,0.0));
1487 0 : put(*vb, -1, true);
1488 0 : break;
1489 0 : case FTMachine::OBSERVED:
1490 : default:
1491 0 : put(*vb, -1, false);
1492 0 : break;
1493 : }
1494 : }
1495 : }
1496 0 : finalizeToSky();
1497 : // Normalize by dividing out weights, etc.
1498 0 : getImage(weight, true);
1499 0 : }
1500 :
1501 4049 : String GridFT::name() const {
1502 :
1503 4049 : return machineName_p;
1504 :
1505 :
1506 : }
1507 :
1508 : }//End of namespace refim
1509 : } //# NAMESPACE CASA - END
1510 :
|
