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