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 22 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType, Float padding,
92 22 : Bool usezero, Bool useDoublePrec)
93 22 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize), tilesize(itilesize),
94 22 : gridder(0), isTiled(false), convType(iconvType),
95 22 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
96 22 : usezero_p(usezero), noPadding_p(false), usePut2_p(false),
97 44 : 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 22 : useDoubleGrid_p=useDoublePrec;
100 : // peek=NULL;
101 22 : }
102 :
103 2534 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
104 2534 : MPosition mLocation, Float padding, Bool usezero, Bool useDoublePrec)
105 2534 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
106 2534 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
107 2534 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
108 5068 : 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 2534 : mLocation_p=mLocation;
111 2534 : tangentSpecified_p=false;
112 2534 : useDoubleGrid_p=useDoublePrec;
113 : // peek=NULL;
114 2534 : }
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 6 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
130 : MPosition mLocation, MDirection mTangent, Float padding,
131 6 : Bool usezero, Bool useDoublePrec)
132 6 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
133 6 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
134 6 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
135 12 : 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 6 : mLocation_p=mLocation;
138 6 : mTangent_p=mTangent;
139 6 : tangentSpecified_p=true;
140 6 : useDoubleGrid_p=useDoublePrec;
141 : // peek=NULL;
142 6 : }
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 267 : GridFT& GridFT::operator=(const GridFT& other)
160 : {
161 267 : if(this!=&other) {
162 : //Do the base parameters
163 267 : FTMachine::operator=(other);
164 :
165 : //private params
166 267 : imageCache=other.imageCache;
167 267 : cachesize=other.cachesize;
168 267 : tilesize=other.tilesize;
169 267 : convType=other.convType;
170 267 : convType.upcase();
171 267 : uvScale.resize();
172 267 : uvOffset.resize();
173 267 : uvScale.assign(other.uvScale);
174 267 : uvOffset.assign(other.uvOffset);
175 267 : if(other.gridder==0)
176 267 : gridder=0;
177 : else{
178 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
179 0 : uvScale, uvOffset,
180 0 : convType);
181 : }
182 267 : convFunc_p.resize();
183 267 : convSupport_p=other.convSupport_p;
184 267 : convSampling_p=other.convSampling_p;
185 267 : convFunc_p=other.convFunc_p;
186 267 : isTiled=other.isTiled;
187 : //lattice=other.lattice;
188 267 : lattice.reset( );
189 267 : tilesize=other.tilesize;
190 267 : arrayLattice.reset( );
191 267 : maxAbsData=other.maxAbsData;
192 267 : centerLoc=other.centerLoc;
193 267 : offsetLoc=other.offsetLoc;
194 267 : padding_p=other.padding_p;
195 267 : usezero_p=other.usezero_p;
196 267 : noPadding_p=other.noPadding_p;
197 267 : machineName_p="GridFT";
198 267 : timemass_p=0.0;
199 267 : timegrid_p=0.0;
200 : // peek = other.peek;
201 : };
202 267 : return *this;
203 : };
204 :
205 : //----------------------------------------------------------------------
206 267 : GridFT::GridFT(const GridFT& other) : FTMachine(), machineName_p("GridFT")
207 : {
208 267 : operator=(other);
209 267 : }
210 :
211 : //-----------------------------------------------------------------------
212 267 : FTMachine* GridFT::cloneFTM(){
213 267 : return new GridFT(*this);
214 : }
215 :
216 : //===============================
217 2197 : Long GridFT::estimateRAM(const CountedPtr<SIImageStore>& imstor){
218 2197 : Long mem=FTMachine::estimateRAM(imstor);
219 2197 : mem=mem*padding_p*padding_p;
220 2197 : return mem;
221 : }
222 :
223 :
224 : //----------------------------------------------------------------------
225 2952 : void GridFT::init() {
226 :
227 2952 : logIO() << LogOrigin("GridFT", "init") << LogIO::NORMAL;
228 :
229 2952 : 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 2952 : isTiled=false;
250 2952 : if(!noPadding_p && padding_p > 1.01){
251 2930 : CompositeNumber cn(uInt(image->shape()(0)*2));
252 2930 : nx = cn.nextLargerEven(Int(padding_p*Float(image->shape()(0))-0.5));
253 2930 : ny = cn.nextLargerEven(Int(padding_p*Float(image->shape()(1))-0.5));
254 2930 : }
255 : else{
256 22 : nx = image->shape()(0);
257 22 : ny = image->shape()(1);
258 : }
259 2952 : npol = image->shape()(2);
260 2952 : nchan = image->shape()(3);
261 : // }
262 :
263 2952 : sumWeight.resize(npol, nchan);
264 :
265 2952 : uvScale.resize(2);
266 2952 : uvScale(0)=(Float(nx)*image->coordinates().increment()(0));
267 2952 : uvScale(1)=(Float(ny)*image->coordinates().increment()(1));
268 2952 : uvOffset.resize(2);
269 2952 : uvOffset(0)=nx/2;
270 2952 : uvOffset(1)=ny/2;
271 :
272 : // Now set up the gridder. The possibilities are BOX and SF
273 2952 : if(gridder) delete gridder; gridder=0;
274 2952 : convType.upcase();
275 5904 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
276 2952 : uvScale, uvOffset,
277 2952 : convType);
278 :
279 2952 : convFunc_p.resize();
280 2952 : convSupport_p=gridder->cSupport()(0);
281 2952 : convSampling_p=gridder->cSampling();
282 2952 : 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 2952 : if(imageCache) delete imageCache; imageCache=0;
288 :
289 2952 : 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 2952 : }
310 :
311 : // This is nasty, we should use CountedPointers here.
312 5658 : GridFT::~GridFT() {
313 2829 : if(imageCache) delete imageCache; imageCache=0;
314 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
315 2829 : if(gridder) delete gridder; gridder=0;
316 5658 : }
317 :
318 : // Initialize for a transform from the Sky domain. This means that
319 : // we grid-correct, and FFT the image
320 :
321 785 : void GridFT::initializeToVis(ImageInterface<Complex>& iimage,
322 : const vi::VisBuffer2& vb)
323 : {
324 785 : image=&iimage;
325 785 : toVis_p=true;
326 :
327 785 : ok();
328 :
329 785 : init();
330 : // peek->reset();
331 : // Initialize the maps for polarization and channel. These maps
332 : // translate visibility indices into image indices
333 785 : 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 785 : prepGridForDegrid();
351 :
352 : //AlwaysAssert(lattice, AipsError);
353 :
354 :
355 :
356 785 : }
357 :
358 785 : void GridFT::prepGridForDegrid(){
359 785 : IPosition gridShape(4, nx, ny, npol, nchan);
360 785 : 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 785 : griddedData.set(Complex(0.0));
366 :
367 785 : IPosition stride(4, 1);
368 1570 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
369 1570 : IPosition trc(blc+image->shape()-stride);
370 :
371 785 : IPosition start(4, 0);
372 785 : griddedData(blc, trc) = image->getSlice(start, image->shape());
373 :
374 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
375 785 : arrayLattice.reset( new ArrayLattice<Complex>(griddedData) );
376 785 : lattice=arrayLattice;
377 : //logIO() << LogIO::DEBUGGING
378 : // << "Starting grid correction and FFT of image" << LogIO::POST;
379 :
380 : // Do the Grid-correction.
381 : {
382 785 : Vector<Complex> correction(nx);
383 785 : correction=Complex(1.0, 0.0);
384 : // Do the Grid-correction
385 785 : IPosition cursorShape(4, nx, 1, 1, 1);
386 785 : IPosition axisPath(4, 0, 1, 2, 3);
387 785 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
388 785 : LatticeIterator<Complex> lix(*lattice, lsx);
389 555009 : for(lix.reset();!lix.atEnd();lix++) {
390 554224 : gridder->correctX1D(correction, lix.position()(1));
391 554224 : lix.rwVectorCursor()/=correction;
392 : }
393 785 : }
394 785 : image->clearCache();
395 : // Now do the FFT2D in place
396 : //LatticeFFT::cfft2d(*lattice);
397 785 : ft_p.c2cFFT(*lattice);
398 : //logIO() << LogIO::DEBUGGING
399 : // << "Finished grid correction and FFT of image" << LogIO::POST;
400 :
401 785 : }
402 :
403 :
404 785 : void GridFT::finalizeToVis()
405 : {
406 :
407 785 : logIO() << LogOrigin("GridFT", "finalizeToVis") << LogIO::NORMAL;
408 785 : logIO() <<LogIO::NORMAL2<< "Time to degrid " << timedegrid_p <<LogIO::POST;
409 785 : timedegrid_p=0.0;
410 :
411 785 : if(arrayLattice) arrayLattice=nullptr;
412 785 : if(lattice) lattice=nullptr;
413 785 : griddedData.resize();
414 785 : 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 785 : }
426 :
427 :
428 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
429 : // grid.
430 2167 : void GridFT::initializeToSky(ImageInterface<Complex>& iimage,
431 : Matrix<Float>& weight, const vi::VisBuffer2& vb)
432 : {
433 : // image always points to the image
434 2167 : image=&iimage;
435 2167 : toVis_p=false;
436 :
437 2167 : init();
438 :
439 : // Initialize the maps for polarization and channel. These maps
440 : // translate visibility indices into image indices
441 2167 : initMaps(vb);
442 :
443 :
444 :
445 2167 : sumWeight=0.0;
446 2167 : weight.resize(sumWeight.shape());
447 2167 : weight=0.0;
448 :
449 :
450 2167 : IPosition gridShape(4, nx, ny, npol, nchan);
451 2167 : if( !useDoubleGrid_p )
452 : {
453 22 : griddedData.resize(gridShape);
454 22 : griddedData=Complex(0.0);
455 : }
456 : else{
457 :
458 2145 : griddedData2.resize(gridShape);
459 2145 : griddedData2=DComplex(0.0);
460 : }
461 2167 : image->clearCache();
462 : //iimage.get(griddedData, false);
463 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
464 :
465 :
466 : //AlwaysAssert(lattice, AipsError);
467 2167 : }
468 :
469 :
470 :
471 2166 : 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 2166 : logIO() << LogOrigin("GridFT", "finalizeToSky") << LogIO::NORMAL;
477 2166 : logIO()<<LogIO::NORMAL2 <<"Time to massage data " << timemass_p << LogIO::POST;
478 2166 : logIO()<< LogIO::NORMAL2 <<"Time to grid data " << timegrid_p << LogIO::POST;
479 2166 : timemass_p=0.0;
480 2166 : timegrid_p=0.0;
481 2166 : 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 2166 : }
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 530188 : void GridFT::put(const vi::VisBuffer2& vb, Int row, Bool dopsf,
684 : FTMachine::Type type)
685 : {
686 :
687 530188 : 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 530187 : 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 530187 : if(max(chanMap)==-1)
707 : {
708 : // peek->reset();
709 136 : return;
710 : }
711 :
712 530051 : Timer tim;
713 530051 : tim.mark();
714 :
715 : //const Matrix<Float> *imagingweight;
716 : //imagingweight=&(vb.imagingWeight());
717 530051 : Matrix<Float> imagingweight;
718 530051 : getImagingWeight(imagingweight, vb);
719 :
720 530051 : if(dopsf) {type=FTMachine::PSF;}
721 :
722 530051 : Cube<Complex> data;
723 : //Fortran gridder need the flag as ints
724 530051 : Cube<Int> flags;
725 530051 : Matrix<Float> elWeight;
726 530051 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
727 :
728 :
729 : Bool iswgtCopy;
730 : const Float *wgtStorage;
731 530051 : wgtStorage=elWeight.getStorage(iswgtCopy);
732 :
733 : Bool isCopy;
734 : const Complex *datStorage;
735 530051 : if(!dopsf)
736 323698 : datStorage=data.getStorage(isCopy);
737 : else
738 206353 : datStorage=0;
739 : // If row is -1 then we pass through all rows
740 : Int startRow, endRow, nRow;
741 530051 : if (row==-1) {
742 530051 : nRow=vb.nRows();
743 530051 : startRow=0;
744 530051 : 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 530051 : const IPosition &fs=flags.shape();
763 530051 : std::vector<Int> s(fs.begin(),fs.end());
764 530051 : Int nvispol=s[0];
765 530051 : Int nvischan=s[1];
766 530051 : Int nvisrow=s[2];
767 530051 : Matrix<Double> uvw(negateUV(vb));
768 :
769 530051 : Vector<Double> dphase(vb.nRows());
770 530051 : Cube<Int> loc(2, nvischan, nRow);
771 530051 : Cube<Int> off(2, nvischan, nRow);
772 530051 : Matrix<Complex> phasor(nvischan, nRow);
773 530051 : Int csamp=convSampling_p;
774 530051 : dphase=0.0;
775 :
776 530051 : timemass_p +=tim.real();
777 530051 : tim.mark();
778 530051 : rotateUVW(uvw, dphase, vb);
779 530051 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
780 : Bool delphase;
781 530051 : Complex * phasorstor=phasor.getStorage(delphase);
782 530051 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
783 530051 : const Double * scalestor=uvScale.getStorage(del);
784 530051 : const Double * offsetstor=uvOffset.getStorage(del);
785 530051 : const Double* uvstor= uvw.getStorage(del);
786 530051 : Int * locstor=loc.getStorage(del);
787 530051 : Int * offstor=off.getStorage(del);
788 530051 : const Double *dpstor=dphase.getStorage(del);
789 : //Vector<Double> f1=interpVisFreq_p.copy();
790 530051 : Int nvchan=nvischan;
791 : Int irow;
792 530051 : Double cinv=Double(1.0)/C::c;
793 530051 : Int dow=0;
794 530051 : Int nth=1;
795 : #ifdef _OPENMP
796 530051 : if(numthreads_p >0){
797 0 : nth=min(numthreads_p, omp_get_max_threads());
798 : }
799 : else{
800 530051 : nth= omp_get_max_threads();
801 : }
802 : //nth=min(4,nth);
803 : #endif
804 :
805 :
806 530051 : #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 530051 : Int idopsf=0;
819 530051 : 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 530051 : Vector<Int> rowFlags(vb.nRows());
831 530051 : rowFlags=0;
832 530051 : rowFlags(vb.flagRow())=true;
833 530051 : if(!usezero_p) {
834 177541900 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
835 177017863 : 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 530051 : Int csupp=convSupport_p;
843 :
844 530051 : const Double * convfuncstor=(convFunc_p).getStorage(del);
845 :
846 : // cerr <<"Poffset " << min(off) << " " << max(off) << " length " << gridder->cFunction().shape() << endl;
847 :
848 :
849 :
850 530051 : ixsub=1;
851 530051 : iysub=1;
852 530051 : if (nth >4){
853 249770 : ixsub=8;
854 249770 : iysub=8;
855 : }
856 280281 : else if(nth >1) {
857 280281 : ixsub=2;
858 280281 : iysub=2;
859 : }
860 :
861 :
862 530051 : Int rbeg=startRow+1;
863 530051 : Int rend=endRow+1;
864 530051 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
865 17636455 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
866 17106404 : sumwgt[icounter].resize(sumWeight.shape());
867 17106404 : sumwgt[icounter].set(0.0);
868 : }
869 530051 : if(doneThreadPartition_p < 0){
870 1200 : xsect_p.resize(ixsub*iysub);
871 1200 : ysect_p.resize(ixsub*iysub);
872 1200 : nxsect_p.resize(ixsub*iysub);
873 1200 : nysect_p.resize(ixsub*iysub);
874 37260 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
875 36060 : findGridSector(nx, ny, ixsub, iysub, 0, 0, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
876 : }
877 : }
878 530051 : Vector<Int> xsect, ysect, nxsect, nysect;
879 530051 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
880 530051 : const Int* pmapstor=polMap.getStorage(del);
881 530051 : const Int *cmapstor=chanMap.getStorage(del);
882 530051 : Int nc=nchan;
883 530051 : Int np=npol;
884 530051 : Int nxp=nx;
885 530051 : Int nyp=ny;
886 530051 : const Int * flagstor=flags.getStorage(del);
887 530051 : const Int * rowflagstor=rowFlags.getStorage(del);
888 : ////////////////////////
889 :
890 : Bool gridcopy;
891 530051 : if(useDoubleGrid_p){
892 524037 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
893 524037 : #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 17577585 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
929 17053548 : sumWeight=sumWeight+sumwgt[icounter];
930 : }
931 : //phasor.putStorage(phasorstor, delphase);
932 524037 : griddedData2.putStorage(gridstor, gridcopy);
933 524037 : if(dopsf && (nth >4))
934 102130 : tweakGridSector(nx, ny, ixsub, iysub);
935 : }
936 : else{
937 6014 : Complex *gridstor=griddedData.getStorage(gridcopy);
938 6014 : #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 58870 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
977 52856 : sumWeight=sumWeight+sumwgt[icounter];
978 : }
979 :
980 6014 : griddedData.putStorage(gridstor, gridcopy);
981 6014 : if(dopsf && (nth > 4))
982 480 : tweakGridSector(nx, ny, ixsub, iysub);
983 : }
984 : // cerr << "sunweight " << sumWeight << endl;
985 :
986 530051 : timegrid_p+=tim.real();
987 :
988 530051 : if(!dopsf)
989 323698 : data.freeStorage(datStorage, isCopy);
990 530051 : elWeight.freeStorage(wgtStorage,iswgtCopy);
991 :
992 : // peek->reset();
993 530051 : }
994 :
995 22 : void GridFT::modifyConvFunc(const Vector<Double>& convFunc, Int convSupport, Int convSampling){
996 22 : convFunc_p.resize();
997 22 : convFunc_p=convFunc;
998 22 : convSupport_p=convSupport;
999 22 : convSampling_p=convSampling;
1000 :
1001 22 : }
1002 :
1003 164612 : void GridFT::get(vi::VisBuffer2& vb, Int row)
1004 : {
1005 :
1006 164612 : gridOk(convSupport_p);
1007 : // If row is -1 then we pass through all rows
1008 : Int startRow, endRow, nRow;
1009 164612 : if (row < 0) {
1010 164612 : nRow=vb.nRows();
1011 164612 : startRow=0;
1012 164612 : 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 164612 : matchChannel(vb);
1026 :
1027 :
1028 : //cerr << "GETchanMap " << chanMap << endl;
1029 : //No point in reading data if its not matching in frequency
1030 164612 : if(max(chanMap)==-1)
1031 58 : return;
1032 :
1033 :
1034 :
1035 : // Get the uvws in a form that Fortran can use
1036 164554 : Matrix<Double> uvw(negateUV(vb));
1037 164554 : Vector<Double> dphase(vb.nRows());
1038 164554 : dphase=0.0;
1039 164554 : rotateUVW(uvw, dphase, vb);
1040 164554 : 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 164554 : Cube<Complex> data;
1053 164554 : Cube<Int> flags;
1054 164554 : 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 164554 : Int nvp=data.shape()(0);
1060 164554 : Int nvc=data.shape()(1);
1061 164554 : Int nvisrow=data.shape()(2);
1062 :
1063 :
1064 : //cerr << "get flags " << min(flags) << " " << max(flags) << endl;
1065 : Complex *datStorage;
1066 : Bool isCopy;
1067 164554 : datStorage=data.getStorage(isCopy);
1068 :
1069 : ///
1070 164554 : Cube<Int> loc(2, nvc, nvisrow);
1071 164554 : Cube<Int> off(2, nvc, nRow);
1072 164554 : Matrix<Complex> phasor(nvc, nRow);
1073 164554 : Int csamp=convSampling_p;
1074 : Bool delphase;
1075 : Bool del;
1076 164554 : Complex * phasorstor=phasor.getStorage(delphase);
1077 164554 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1078 164554 : const Double * scalestor=uvScale.getStorage(del);
1079 164554 : const Double * offsetstor=uvOffset.getStorage(del);
1080 164554 : const Double* uvstor= uvw.getStorage(del);
1081 164554 : Int * locstor=loc.getStorage(del);
1082 164554 : Int * offstor=off.getStorage(del);
1083 164554 : const Double *dpstor=dphase.getStorage(del);
1084 : //Vector<Double> f1=interpVisFreq_p.copy();
1085 164554 : Int nvchan=nvc;
1086 : Int irow;
1087 164554 : Double cinv=Double(1.0)/C::c;
1088 164554 : Int dow=0;
1089 164554 : Int nth=1;
1090 : #ifdef _OPENMP
1091 164554 : if(numthreads_p >0){
1092 0 : nth=min(numthreads_p, omp_get_max_threads());
1093 : }
1094 : else{
1095 164554 : nth= omp_get_max_threads();
1096 : }
1097 : //nth=min(4,nth);
1098 : #endif
1099 :
1100 :
1101 164554 : #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 164554 : Int rbeg=startRow+1;
1116 164554 : Int rend=endRow+1;
1117 :
1118 :
1119 164554 : Vector<Int> rowFlags(vb.nRows());
1120 164554 : rowFlags=0;
1121 164554 : rowFlags(vb.flagRow())=true;
1122 : //cerr << "rowFlags " << rowFlags << endl;
1123 164554 : if(!usezero_p) {
1124 56762618 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1125 56598064 : 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 164554 : const Complex* gridstor=griddedData.getStorage(delgrid);
1135 164554 : const Double * convfuncstor=(convFunc_p).getStorage(del);
1136 :
1137 164554 : const Int* pmapstor=polMap.getStorage(del);
1138 164554 : const Int *cmapstor=chanMap.getStorage(del);
1139 164554 : Int nc=nchan;
1140 164554 : Int np=npol;
1141 164554 : Int nxp=nx;
1142 164554 : Int nyp=ny;
1143 164554 : Int csupp=convSupport_p;
1144 164554 : const Int * flagstor=flags.getStorage(del);
1145 164554 : const Int * rowflagstor=rowFlags.getStorage(del);
1146 :
1147 :
1148 164554 : Int npart=nth;
1149 :
1150 :
1151 164554 : Int ix=0;
1152 164554 : #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 164554 : data.putStorage(datStorage, isCopy);
1180 164554 : griddedData.freeStorage(gridstor, delgrid);
1181 : //cerr << "Get min max " << min(data) << " " << max(data) << endl;
1182 : }
1183 164554 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1184 :
1185 164554 : }
1186 :
1187 :
1188 :
1189 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1190 : // return the resulting image
1191 2144 : ImageInterface<Complex>& GridFT::getImage(Matrix<Float>& weights, Bool normalize)
1192 : {
1193 : //AlwaysAssert(lattice, AipsError);
1194 2144 : AlwaysAssert(gridder, AipsError);
1195 2144 : AlwaysAssert(image, AipsError);
1196 2144 : logIO() << LogOrigin("GridFT", "getImage") << LogIO::NORMAL;
1197 :
1198 2144 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
1199 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
1200 : }
1201 2144 : weights.resize(sumWeight.shape());
1202 :
1203 2144 : convertArray(weights, sumWeight);
1204 : // If the weights are all zero then we cannot normalize
1205 : // otherwise we don't care.
1206 2144 : 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 2144 : if(useDoubleGrid_p)
1233 : {
1234 2144 : ArrayLattice<DComplex> darrayLattice(griddedData2);
1235 : //LatticeFFT::cfft2d(darrayLattice,false);
1236 2144 : ft_p.c2cFFT(darrayLattice, False);
1237 2144 : griddedData.resize(griddedData2.shape());
1238 2144 : convertArray(griddedData, griddedData2);
1239 :
1240 2144 : SynthesisUtilMethods::getResource("mem peak in getImage");
1241 :
1242 : //Don't need the double-prec grid anymore...
1243 2144 : griddedData2.resize();
1244 2144 : arrayLattice.reset( new ArrayLattice<Complex>(griddedData) );
1245 2144 : lattice=arrayLattice;
1246 2144 : }
1247 : else{
1248 0 : arrayLattice.reset( new ArrayLattice<Complex>(griddedData) );
1249 0 : lattice=arrayLattice;
1250 :
1251 : //LatticeFFT::cfft2d(*lattice,false);
1252 0 : ft_p.c2cFFT(*lattice, False);
1253 : }
1254 :
1255 :
1256 :
1257 : {
1258 2144 : Int inx = lattice->shape()(0);
1259 2144 : Int iny = lattice->shape()(1);
1260 2144 : Vector<Complex> correction(inx);
1261 2144 : correction=Complex(1.0, 0.0);
1262 : // Do the Grid-correction
1263 2144 : IPosition cursorShape(4, inx, 1, 1, 1);
1264 2144 : IPosition axisPath(4, 0, 1, 2, 3);
1265 2144 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1266 2144 : LatticeIterator<Complex> lix(*lattice, lsx);
1267 1583072 : for(lix.reset();!lix.atEnd();lix++) {
1268 1580928 : Int pol=lix.position()(2);
1269 1580928 : Int chan=lix.position()(3);
1270 1580928 : if(weights(pol, chan)!=0.0) {
1271 1523400 : gridder->correctX1D(correction, lix.position()(1));
1272 1523400 : lix.rwVectorCursor()/=correction;
1273 1523400 : if(normalize) {
1274 0 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
1275 0 : lix.rwCursor()*=rnorm;
1276 : }
1277 : else {
1278 1523400 : Complex rnorm(Float(inx)*Float(iny));
1279 1523400 : lix.rwCursor()*=rnorm;
1280 : }
1281 : }
1282 : else {
1283 57528 : lix.woCursor()=0.0;
1284 : }
1285 : }
1286 2144 : }
1287 2144 : LatticeLocker lock1 (*(image), FileLocker::Write);
1288 2144 : if(!isTiled) {
1289 : // Check the section from the image BEFORE converting to a lattice
1290 4288 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1291 2144 : IPosition stride(4, 1);
1292 4288 : IPosition trc(blc+image->shape()-stride);
1293 : // Do the copy
1294 2144 : IPosition start(4, 0);
1295 :
1296 2144 : image->put(griddedData(blc, trc));
1297 2144 : }
1298 2144 : }
1299 2144 : image->flush();
1300 2144 : image->clearCache();
1301 :
1302 2144 : if(arrayLattice) arrayLattice=nullptr;
1303 2144 : if(lattice) lattice=nullptr;
1304 2144 : griddedData.resize();
1305 2144 : 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 45 : 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 45 : Bool retval = true;
1342 45 : Float elpadd=padding_p;
1343 45 : if(toVis_p && withImage)
1344 19 : elpadd=1.0;
1345 : //save the base class variables
1346 45 : 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 45 : outRec.define("cachesize", Int64(cachesize));
1353 45 : outRec.define("tilesize", tilesize);
1354 45 : outRec.define("convtype", convType);
1355 45 : outRec.define("uvscale", uvScale);
1356 45 : outRec.define("uvoffset", uvOffset);
1357 45 : outRec.define("istiled", isTiled);
1358 45 : outRec.define("maxabsdata", maxAbsData);
1359 45 : Vector<Int> center_loc(4), offset_loc(4);
1360 225 : for (Int k=0; k<4 ; k++){
1361 180 : center_loc(k)=centerLoc(k);
1362 180 : offset_loc(k)=offsetLoc(k);
1363 : }
1364 45 : outRec.define("centerloc", center_loc);
1365 45 : outRec.define("offsetloc", offset_loc);
1366 45 : outRec.define("padding", elpadd);
1367 45 : outRec.define("usezero", usezero_p);
1368 45 : outRec.define("nopadding", noPadding_p);
1369 45 : return retval;
1370 45 : }
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 1392947 : void GridFT::ok() {
1431 1392947 : AlwaysAssert(image, AipsError);
1432 1392947 : }
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 4148 : String GridFT::name() const {
1502 :
1503 4148 : return machineName_p;
1504 :
1505 :
1506 : }
1507 :
1508 : }//End of namespace refim
1509 : } //# NAMESPACE CASA - END
1510 :
|
