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: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <msvis/MSVis/VisibilityIterator2.h>
29 : #include <casacore/casa/Quanta/UnitMap.h>
30 : #include <casacore/casa/Quanta/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 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
87 : 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 14 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType, Float padding,
92 14 : Bool usezero, Bool useDoublePrec)
93 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize), tilesize(itilesize),
94 : gridder(0), isTiled(false), convType(iconvType),
95 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
96 : usezero_p(usezero), noPadding_p(false), usePut2_p(false),
97 14 : 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 14 : useDoubleGrid_p=useDoublePrec;
100 : // peek=NULL;
101 14 : }
102 :
103 2432 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
104 2432 : MPosition mLocation, Float padding, Bool usezero, Bool useDoublePrec)
105 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
106 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
107 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
108 2432 : 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 2432 : mLocation_p=mLocation;
111 2432 : tangentSpecified_p=false;
112 2432 : useDoubleGrid_p=useDoublePrec;
113 : // peek=NULL;
114 2432 : }
115 :
116 0 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
117 0 : MDirection mTangent, Float padding, Bool usezero, Bool useDoublePrec)
118 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
119 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
120 : 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 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
133 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
134 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
135 6 : 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 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 2111 : Long GridFT::estimateRAM(const CountedPtr<SIImageStore>& imstor){
218 2111 : Long mem=FTMachine::estimateRAM(imstor);
219 2111 : mem=mem*padding_p*padding_p;
220 2111 : return mem;
221 : }
222 :
223 :
224 : //----------------------------------------------------------------------
225 2842 : void GridFT::init() {
226 :
227 2842 : logIO() << LogOrigin("GridFT", "init") << LogIO::NORMAL;
228 :
229 2842 : 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 2842 : isTiled=false;
250 2842 : if(!noPadding_p && padding_p > 1.01){
251 2828 : CompositeNumber cn(uInt(image->shape()(0)*2));
252 2828 : nx = cn.nextLargerEven(Int(padding_p*Float(image->shape()(0))-0.5));
253 2828 : ny = cn.nextLargerEven(Int(padding_p*Float(image->shape()(1))-0.5));
254 : }
255 : else{
256 14 : nx = image->shape()(0);
257 14 : ny = image->shape()(1);
258 : }
259 2842 : npol = image->shape()(2);
260 2842 : nchan = image->shape()(3);
261 : // }
262 :
263 2842 : sumWeight.resize(npol, nchan);
264 :
265 2842 : uvScale.resize(2);
266 2842 : uvScale(0)=(Float(nx)*image->coordinates().increment()(0));
267 2842 : uvScale(1)=(Float(ny)*image->coordinates().increment()(1));
268 2842 : uvOffset.resize(2);
269 2842 : uvOffset(0)=nx/2;
270 2842 : uvOffset(1)=ny/2;
271 :
272 : // Now set up the gridder. The possibilities are BOX and SF
273 2842 : if(gridder) delete gridder; gridder=0;
274 2842 : convType.upcase();
275 5684 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
276 2842 : uvScale, uvOffset,
277 2842 : convType);
278 :
279 2842 : convFunc_p.resize();
280 2842 : convSupport_p=gridder->cSupport()(0);
281 2842 : convSampling_p=gridder->cSampling();
282 2842 : 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 2842 : if(imageCache) delete imageCache; imageCache=0;
288 :
289 2842 : 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 : }
309 2842 : }
310 :
311 : // This is nasty, we should use CountedPointers here.
312 5432 : GridFT::~GridFT() {
313 2716 : if(imageCache) delete imageCache; imageCache=0;
314 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
315 2716 : if(gridder) delete gridder; gridder=0;
316 5432 : }
317 :
318 : // Initialize for a transform from the Sky domain. This means that
319 : // we grid-correct, and FFT the image
320 :
321 764 : void GridFT::initializeToVis(ImageInterface<Complex>& iimage,
322 : const vi::VisBuffer2& vb)
323 : {
324 764 : image=&iimage;
325 764 : toVis_p=true;
326 :
327 764 : ok();
328 :
329 764 : init();
330 : // peek->reset();
331 : // Initialize the maps for polarization and channel. These maps
332 : // translate visibility indices into image indices
333 764 : 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 764 : prepGridForDegrid();
351 :
352 : //AlwaysAssert(lattice, AipsError);
353 :
354 :
355 :
356 764 : }
357 :
358 764 : void GridFT::prepGridForDegrid(){
359 1528 : IPosition gridShape(4, nx, ny, npol, nchan);
360 764 : 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 764 : griddedData.set(Complex(0.0));
366 :
367 1528 : IPosition stride(4, 1);
368 2292 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
369 2292 : IPosition trc(blc+image->shape()-stride);
370 :
371 1528 : IPosition start(4, 0);
372 764 : griddedData(blc, trc) = image->getSlice(start, image->shape());
373 :
374 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
375 764 : arrayLattice.reset( new ArrayLattice<Complex>(griddedData) );
376 764 : lattice=arrayLattice;
377 : //logIO() << LogIO::DEBUGGING
378 : // << "Starting grid correction and FFT of image" << LogIO::POST;
379 :
380 : // Do the Grid-correction.
381 : {
382 1528 : Vector<Complex> correction(nx);
383 764 : correction=Complex(1.0, 0.0);
384 : // Do the Grid-correction
385 1528 : IPosition cursorShape(4, nx, 1, 1, 1);
386 1528 : IPosition axisPath(4, 0, 1, 2, 3);
387 1528 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
388 1528 : LatticeIterator<Complex> lix(*lattice, lsx);
389 531028 : for(lix.reset();!lix.atEnd();lix++) {
390 530264 : gridder->correctX1D(correction, lix.position()(1));
391 530264 : lix.rwVectorCursor()/=correction;
392 : }
393 : }
394 764 : image->clearCache();
395 : // Now do the FFT2D in place
396 : //LatticeFFT::cfft2d(*lattice);
397 764 : ft_p.c2cFFT(*lattice);
398 : //logIO() << LogIO::DEBUGGING
399 : // << "Finished grid correction and FFT of image" << LogIO::POST;
400 :
401 764 : }
402 :
403 :
404 764 : void GridFT::finalizeToVis()
405 : {
406 :
407 764 : logIO() << LogOrigin("GridFT", "finalizeToVis") << LogIO::NORMAL;
408 764 : logIO() <<LogIO::NORMAL2<< "Time to degrid " << timedegrid_p <<LogIO::POST;
409 764 : timedegrid_p=0.0;
410 :
411 764 : if(arrayLattice) arrayLattice=nullptr;
412 764 : if(lattice) lattice=nullptr;
413 764 : griddedData.resize();
414 764 : 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 : }
425 764 : }
426 :
427 :
428 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
429 : // grid.
430 2078 : void GridFT::initializeToSky(ImageInterface<Complex>& iimage,
431 : Matrix<Float>& weight, const vi::VisBuffer2& vb)
432 : {
433 : // image always points to the image
434 2078 : image=&iimage;
435 2078 : toVis_p=false;
436 :
437 2078 : init();
438 :
439 : // Initialize the maps for polarization and channel. These maps
440 : // translate visibility indices into image indices
441 2078 : initMaps(vb);
442 :
443 :
444 :
445 2078 : sumWeight=0.0;
446 2078 : weight.resize(sumWeight.shape());
447 2078 : weight=0.0;
448 :
449 :
450 4156 : IPosition gridShape(4, nx, ny, npol, nchan);
451 2078 : if( !useDoubleGrid_p )
452 : {
453 14 : griddedData.resize(gridShape);
454 14 : griddedData=Complex(0.0);
455 : }
456 : else{
457 :
458 2064 : griddedData2.resize(gridShape);
459 2064 : griddedData2=DComplex(0.0);
460 : }
461 2078 : image->clearCache();
462 : //iimage.get(griddedData, false);
463 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
464 :
465 :
466 : //AlwaysAssert(lattice, AipsError);
467 2078 : }
468 :
469 :
470 :
471 2077 : 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 2077 : logIO() << LogOrigin("GridFT", "finalizeToSky") << LogIO::NORMAL;
477 2077 : logIO()<<LogIO::NORMAL2 <<"Time to massage data " << timemass_p << LogIO::POST;
478 2077 : logIO()<< LogIO::NORMAL2 <<"Time to grid data " << timegrid_p << LogIO::POST;
479 2077 : timemass_p=0.0;
480 2077 : timegrid_p=0.0;
481 2077 : 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 : }
491 : // peek->terminate();
492 : // peek->reset();
493 2077 : }
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 501698 : void GridFT::put(const vi::VisBuffer2& vb, Int row, Bool dopsf,
684 : FTMachine::Type type)
685 : {
686 :
687 501698 : 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 501697 : 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 501697 : if(max(chanMap)==-1)
707 : {
708 : // peek->reset();
709 136 : return;
710 : }
711 :
712 501561 : Timer tim;
713 501561 : tim.mark();
714 :
715 : //const Matrix<Float> *imagingweight;
716 : //imagingweight=&(vb.imagingWeight());
717 1003122 : Matrix<Float> imagingweight;
718 501561 : getImagingWeight(imagingweight, vb);
719 :
720 501561 : if(dopsf) {type=FTMachine::PSF;}
721 :
722 1003122 : Cube<Complex> data;
723 : //Fortran gridder need the flag as ints
724 1003122 : Cube<Int> flags;
725 1003122 : Matrix<Float> elWeight;
726 501561 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
727 :
728 :
729 : Bool iswgtCopy;
730 : const Float *wgtStorage;
731 501561 : wgtStorage=elWeight.getStorage(iswgtCopy);
732 :
733 : Bool isCopy;
734 : const Complex *datStorage;
735 501561 : if(!dopsf)
736 306684 : datStorage=data.getStorage(isCopy);
737 : else
738 194877 : datStorage=0;
739 : // If row is -1 then we pass through all rows
740 : Int startRow, endRow, nRow;
741 501561 : if (row==-1) {
742 501561 : nRow=vb.nRows();
743 501561 : startRow=0;
744 501561 : 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 501561 : const IPosition &fs=flags.shape();
763 1003122 : std::vector<Int> s(fs.begin(),fs.end());
764 501561 : Int nvispol=s[0];
765 501561 : Int nvischan=s[1];
766 501561 : Int nvisrow=s[2];
767 1003122 : Matrix<Double> uvw(negateUV(vb));
768 :
769 1003122 : Vector<Double> dphase(vb.nRows());
770 1003122 : Cube<Int> loc(2, nvischan, nRow);
771 1003122 : Cube<Int> off(2, nvischan, nRow);
772 1003122 : Matrix<Complex> phasor(nvischan, nRow);
773 501561 : Int csamp=convSampling_p;
774 501561 : dphase=0.0;
775 :
776 501561 : timemass_p +=tim.real();
777 501561 : tim.mark();
778 501561 : rotateUVW(uvw, dphase, vb);
779 501561 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
780 : Bool delphase;
781 501561 : Complex * phasorstor=phasor.getStorage(delphase);
782 501561 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
783 501561 : const Double * scalestor=uvScale.getStorage(del);
784 501561 : const Double * offsetstor=uvOffset.getStorage(del);
785 501561 : const Double* uvstor= uvw.getStorage(del);
786 501561 : Int * locstor=loc.getStorage(del);
787 501561 : Int * offstor=off.getStorage(del);
788 501561 : const Double *dpstor=dphase.getStorage(del);
789 : //Vector<Double> f1=interpVisFreq_p.copy();
790 501561 : Int nvchan=nvischan;
791 : Int irow;
792 501561 : Double cinv=Double(1.0)/C::c;
793 501561 : Int dow=0;
794 501561 : Int nth=1;
795 : #ifdef _OPENMP
796 501561 : if(numthreads_p >0){
797 0 : nth=min(numthreads_p, omp_get_max_threads());
798 : }
799 : else{
800 501561 : nth= omp_get_max_threads();
801 : }
802 : //nth=min(4,nth);
803 : #endif
804 :
805 :
806 501561 : #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 501561 : Int idopsf=0;
819 501561 : 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 1003122 : Vector<Int> rowFlags(vb.nRows());
831 501561 : rowFlags=0;
832 501561 : rowFlags(vb.flagRow())=true;
833 501561 : if(!usezero_p) {
834 168987676 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
835 168491509 : 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 501561 : Int csupp=convSupport_p;
843 :
844 501561 : const Double * convfuncstor=(convFunc_p).getStorage(del);
845 :
846 : // cerr <<"Poffset " << min(off) << " " << max(off) << " length " << gridder->cFunction().shape() << endl;
847 :
848 :
849 :
850 501561 : ixsub=1;
851 501561 : iysub=1;
852 501561 : if (nth >4){
853 0 : ixsub=8;
854 0 : iysub=8;
855 : }
856 501561 : else if(nth >1) {
857 0 : ixsub=2;
858 0 : iysub=2;
859 : }
860 :
861 :
862 501561 : Int rbeg=startRow+1;
863 501561 : Int rend=endRow+1;
864 1003122 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
865 1003122 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
866 501561 : sumwgt[icounter].resize(sumWeight.shape());
867 501561 : sumwgt[icounter].set(0.0);
868 : }
869 501561 : if(doneThreadPartition_p < 0){
870 1146 : xsect_p.resize(ixsub*iysub);
871 1146 : ysect_p.resize(ixsub*iysub);
872 1146 : nxsect_p.resize(ixsub*iysub);
873 1146 : nysect_p.resize(ixsub*iysub);
874 2292 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
875 1146 : findGridSector(nx, ny, ixsub, iysub, 0, 0, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
876 : }
877 : }
878 1003122 : Vector<Int> xsect, ysect, nxsect, nysect;
879 501561 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
880 501561 : const Int* pmapstor=polMap.getStorage(del);
881 501561 : const Int *cmapstor=chanMap.getStorage(del);
882 501561 : Int nc=nchan;
883 501561 : Int np=npol;
884 501561 : Int nxp=nx;
885 501561 : Int nyp=ny;
886 501561 : const Int * flagstor=flags.getStorage(del);
887 501561 : const Int * rowflagstor=rowFlags.getStorage(del);
888 : ////////////////////////
889 :
890 : Bool gridcopy;
891 501561 : if(useDoubleGrid_p){
892 496167 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
893 496167 : #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 992334 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
929 496167 : sumWeight=sumWeight+sumwgt[icounter];
930 : }
931 : //phasor.putStorage(phasorstor, delphase);
932 496167 : griddedData2.putStorage(gridstor, gridcopy);
933 496167 : if(dopsf && (nth >4))
934 0 : tweakGridSector(nx, ny, ixsub, iysub);
935 : }
936 : else{
937 5394 : Complex *gridstor=griddedData.getStorage(gridcopy);
938 5394 : #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 10788 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
977 5394 : sumWeight=sumWeight+sumwgt[icounter];
978 : }
979 :
980 5394 : griddedData.putStorage(gridstor, gridcopy);
981 5394 : if(dopsf && (nth > 4))
982 0 : tweakGridSector(nx, ny, ixsub, iysub);
983 : }
984 : // cerr << "sunweight " << sumWeight << endl;
985 :
986 501561 : timegrid_p+=tim.real();
987 :
988 501561 : if(!dopsf)
989 306684 : data.freeStorage(datStorage, isCopy);
990 501561 : elWeight.freeStorage(wgtStorage,iswgtCopy);
991 :
992 : // peek->reset();
993 : }
994 :
995 14 : void GridFT::modifyConvFunc(const Vector<Double>& convFunc, Int convSupport, Int convSampling){
996 14 : convFunc_p.resize();
997 14 : convFunc_p=convFunc;
998 14 : convSupport_p=convSupport;
999 14 : convSampling_p=convSampling;
1000 :
1001 14 : }
1002 :
1003 157782 : void GridFT::get(vi::VisBuffer2& vb, Int row)
1004 : {
1005 :
1006 157782 : gridOk(convSupport_p);
1007 : // If row is -1 then we pass through all rows
1008 : Int startRow, endRow, nRow;
1009 157782 : if (row < 0) {
1010 157782 : nRow=vb.nRows();
1011 157782 : startRow=0;
1012 157782 : 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 157782 : matchChannel(vb);
1026 :
1027 :
1028 : //cerr << "GETchanMap " << chanMap << endl;
1029 : //No point in reading data if its not matching in frequency
1030 157782 : if(max(chanMap)==-1)
1031 58 : return;
1032 :
1033 :
1034 :
1035 : // Get the uvws in a form that Fortran can use
1036 315448 : Matrix<Double> uvw(negateUV(vb));
1037 315448 : Vector<Double> dphase(vb.nRows());
1038 157724 : dphase=0.0;
1039 157724 : rotateUVW(uvw, dphase, vb);
1040 157724 : 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 315448 : Cube<Complex> data;
1053 315448 : Cube<Int> flags;
1054 157724 : 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 157724 : Int nvp=data.shape()(0);
1060 157724 : Int nvc=data.shape()(1);
1061 157724 : Int nvisrow=data.shape()(2);
1062 :
1063 :
1064 : //cerr << "get flags " << min(flags) << " " << max(flags) << endl;
1065 : Complex *datStorage;
1066 : Bool isCopy;
1067 157724 : datStorage=data.getStorage(isCopy);
1068 :
1069 : ///
1070 315448 : Cube<Int> loc(2, nvc, nvisrow);
1071 315448 : Cube<Int> off(2, nvc, nRow);
1072 315448 : Matrix<Complex> phasor(nvc, nRow);
1073 157724 : Int csamp=convSampling_p;
1074 : Bool delphase;
1075 : Bool del;
1076 157724 : Complex * phasorstor=phasor.getStorage(delphase);
1077 157724 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1078 157724 : const Double * scalestor=uvScale.getStorage(del);
1079 157724 : const Double * offsetstor=uvOffset.getStorage(del);
1080 157724 : const Double* uvstor= uvw.getStorage(del);
1081 157724 : Int * locstor=loc.getStorage(del);
1082 157724 : Int * offstor=off.getStorage(del);
1083 157724 : const Double *dpstor=dphase.getStorage(del);
1084 : //Vector<Double> f1=interpVisFreq_p.copy();
1085 157724 : Int nvchan=nvc;
1086 : Int irow;
1087 157724 : Double cinv=Double(1.0)/C::c;
1088 157724 : Int dow=0;
1089 157724 : Int nth=1;
1090 : #ifdef _OPENMP
1091 157724 : if(numthreads_p >0){
1092 0 : nth=min(numthreads_p, omp_get_max_threads());
1093 : }
1094 : else{
1095 157724 : nth= omp_get_max_threads();
1096 : }
1097 : //nth=min(4,nth);
1098 : #endif
1099 :
1100 :
1101 157724 : #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 157724 : Int rbeg=startRow+1;
1116 157724 : Int rend=endRow+1;
1117 :
1118 :
1119 315448 : Vector<Int> rowFlags(vb.nRows());
1120 157724 : rowFlags=0;
1121 157724 : rowFlags(vb.flagRow())=true;
1122 : //cerr << "rowFlags " << rowFlags << endl;
1123 157724 : if(!usezero_p) {
1124 54180610 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1125 54022886 : 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 157724 : const Complex* gridstor=griddedData.getStorage(delgrid);
1135 157724 : const Double * convfuncstor=(convFunc_p).getStorage(del);
1136 :
1137 157724 : const Int* pmapstor=polMap.getStorage(del);
1138 157724 : const Int *cmapstor=chanMap.getStorage(del);
1139 157724 : Int nc=nchan;
1140 157724 : Int np=npol;
1141 157724 : Int nxp=nx;
1142 157724 : Int nyp=ny;
1143 157724 : Int csupp=convSupport_p;
1144 157724 : const Int * flagstor=flags.getStorage(del);
1145 157724 : const Int * rowflagstor=rowFlags.getStorage(del);
1146 :
1147 :
1148 157724 : Int npart=nth;
1149 :
1150 :
1151 157724 : Int ix=0;
1152 157724 : #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 157724 : data.putStorage(datStorage, isCopy);
1180 157724 : griddedData.freeStorage(gridstor, delgrid);
1181 : //cerr << "Get min max " << min(data) << " " << max(data) << endl;
1182 : }
1183 157724 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1184 :
1185 : }
1186 :
1187 :
1188 :
1189 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1190 : // return the resulting image
1191 2063 : ImageInterface<Complex>& GridFT::getImage(Matrix<Float>& weights, Bool normalize)
1192 : {
1193 : //AlwaysAssert(lattice, AipsError);
1194 2063 : AlwaysAssert(gridder, AipsError);
1195 2063 : AlwaysAssert(image, AipsError);
1196 2063 : logIO() << LogOrigin("GridFT", "getImage") << LogIO::NORMAL;
1197 :
1198 2063 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
1199 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
1200 : }
1201 2063 : weights.resize(sumWeight.shape());
1202 :
1203 2063 : convertArray(weights, sumWeight);
1204 : // If the weights are all zero then we cannot normalize
1205 : // otherwise we don't care.
1206 2063 : 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 2063 : if(useDoubleGrid_p)
1233 : {
1234 4126 : ArrayLattice<DComplex> darrayLattice(griddedData2);
1235 : //LatticeFFT::cfft2d(darrayLattice,false);
1236 2063 : ft_p.c2cFFT(darrayLattice, False);
1237 2063 : griddedData.resize(griddedData2.shape());
1238 2063 : convertArray(griddedData, griddedData2);
1239 :
1240 2063 : SynthesisUtilMethods::getResource("mem peak in getImage");
1241 :
1242 : //Don't need the double-prec grid anymore...
1243 2063 : griddedData2.resize();
1244 2063 : arrayLattice.reset( new ArrayLattice<Complex>(griddedData) );
1245 2063 : lattice=arrayLattice;
1246 : }
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 2063 : Int inx = lattice->shape()(0);
1259 2063 : Int iny = lattice->shape()(1);
1260 4126 : Vector<Complex> correction(inx);
1261 2063 : correction=Complex(1.0, 0.0);
1262 : // Do the Grid-correction
1263 4126 : IPosition cursorShape(4, inx, 1, 1, 1);
1264 4126 : IPosition axisPath(4, 0, 1, 2, 3);
1265 4126 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1266 4126 : LatticeIterator<Complex> lix(*lattice, lsx);
1267 1510951 : for(lix.reset();!lix.atEnd();lix++) {
1268 1508888 : Int pol=lix.position()(2);
1269 1508888 : Int chan=lix.position()(3);
1270 1508888 : if(weights(pol, chan)!=0.0) {
1271 1462880 : gridder->correctX1D(correction, lix.position()(1));
1272 1462880 : lix.rwVectorCursor()/=correction;
1273 1462880 : if(normalize) {
1274 0 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
1275 0 : lix.rwCursor()*=rnorm;
1276 : }
1277 : else {
1278 1462880 : Complex rnorm(Float(inx)*Float(iny));
1279 1462880 : lix.rwCursor()*=rnorm;
1280 : }
1281 : }
1282 : else {
1283 46008 : lix.woCursor()=0.0;
1284 : }
1285 : }
1286 : }
1287 4126 : LatticeLocker lock1 (*(image), FileLocker::Write);
1288 2063 : if(!isTiled) {
1289 : // Check the section from the image BEFORE converting to a lattice
1290 6189 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1291 4126 : IPosition stride(4, 1);
1292 6189 : IPosition trc(blc+image->shape()-stride);
1293 : // Do the copy
1294 2063 : IPosition start(4, 0);
1295 :
1296 2063 : image->put(griddedData(blc, trc));
1297 : }
1298 : }
1299 2063 : image->flush();
1300 2063 : image->clearCache();
1301 :
1302 2063 : if(arrayLattice) arrayLattice=nullptr;
1303 2063 : if(lattice) lattice=nullptr;
1304 2063 : griddedData.resize();
1305 2063 : 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 84 : 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 : }
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 : }
1429 :
1430 1322176 : void GridFT::ok() {
1431 1322176 : AlwaysAssert(image, AipsError);
1432 1322176 : }
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 3975 : String GridFT::name() const {
1502 :
1503 3975 : return machineName_p;
1504 :
1505 :
1506 : }
1507 :
1508 : }//End of namespace refim
1509 : } //# NAMESPACE CASA - END
1510 :
|
