Line data Source code
1 : //# GridFT.cc: Implementation of GridFT class
2 : //# Copyright (C) 1997-2016
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <msvis/MSVis/VisibilityIterator2.h>
29 : #include <casacore/casa/Quanta/UnitMap.h>
30 : #include <casacore/casa/Quanta/UnitVal.h>
31 : #include <casacore/measures/Measures/Stokes.h>
32 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
33 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
34 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
35 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
36 : #include <casacore/coordinates/Coordinates/Projection.h>
37 : #include <casacore/ms/MeasurementSets/MSColumns.h>
38 : #include <casacore/casa/BasicSL/Constants.h>
39 : #include <casacore/scimath/Mathematics/FFTServer.h>
40 : #include <casacore/scimath/Mathematics/RigidVector.h>
41 : #include <msvis/MSVis/StokesVector.h>
42 : #include <synthesis/TransformMachines/StokesImageUtil.h>
43 : #include <synthesis/TransformMachines2/GridFT.h>
44 : #include <synthesis/Utilities/FFT2D.h>
45 : #include <msvis/MSVis/VisBuffer2.h>
46 : #include <casacore/images/Images/ImageInterface.h>
47 : #include <casacore/images/Images/PagedImage.h>
48 : #include <casacore/casa/Containers/Block.h>
49 : #include <casacore/casa/Containers/Record.h>
50 : #include <casacore/casa/Arrays/ArrayLogical.h>
51 : #include <casacore/casa/Arrays/ArrayMath.h>
52 : #include <casacore/casa/Arrays/Array.h>
53 : #include <casacore/casa/Arrays/MaskedArray.h>
54 : #include <casacore/casa/Arrays/Vector.h>
55 : #include <casacore/casa/Arrays/Matrix.h>
56 : #include <casacore/casa/Arrays/Cube.h>
57 : #include <casacore/casa/Arrays/MatrixIter.h>
58 : #include <casacore/casa/BasicSL/String.h>
59 : #include <casacore/casa/Utilities/Assert.h>
60 : #include <casacore/casa/Exceptions/Error.h>
61 : #include <casacore/lattices/Lattices/ArrayLattice.h>
62 : #include <casacore/measures/Measures/UVWMachine.h>
63 : #include <casacore/lattices/Lattices/SubLattice.h>
64 : #include <casacore/lattices/LRegions/LCBox.h>
65 : #include <casacore/lattices/Lattices/LatticeCache.h>
66 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
67 : #include <casacore/lattices/Lattices/LatticeIterator.h>
68 : #include <casacore/lattices/Lattices/LatticeStepper.h>
69 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
70 : #include <casacore/casa/Utilities/CompositeNumber.h>
71 : #include <casacore/casa/OS/Timer.h>
72 : #include <sstream>
73 : #ifdef _OPENMP
74 : #include <omp.h>
75 : #endif
76 :
77 : using namespace casacore;
78 : namespace casa { //# NAMESPACE CASA - BEGIN
79 : namespace refim {//# namespace for imaging refactor
80 : using namespace casacore;
81 : using namespace casa;
82 : using namespace casacore;
83 : using namespace casa::refim;
84 :
85 0 : GridFT::GridFT() : FTMachine(), padding_p(1.0), imageCache(0), cachesize(1000000), tilesize(1000), gridder(0), isTiled(false), convType("SF"),
86 0 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
87 0 : usezero_p(false), noPadding_p(false), usePut2_p(false),
88 0 : machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0){
89 :
90 0 : }
91 22 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType, Float padding,
92 22 : Bool usezero, Bool useDoublePrec)
93 22 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize), tilesize(itilesize),
94 22 : gridder(0), isTiled(false), convType(iconvType),
95 22 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
96 22 : usezero_p(usezero), noPadding_p(false), usePut2_p(false),
97 44 : machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
98 : {
99 22 : useDoubleGrid_p=useDoublePrec;
100 : // peek=NULL;
101 22 : }
102 :
103 2530 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
104 2530 : MPosition mLocation, Float padding, Bool usezero, Bool useDoublePrec)
105 2530 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
106 2530 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
107 2530 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
108 5060 : 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 2530 : mLocation_p=mLocation;
111 2530 : tangentSpecified_p=false;
112 2530 : useDoubleGrid_p=useDoublePrec;
113 : // peek=NULL;
114 2530 : }
115 :
116 0 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
117 0 : MDirection mTangent, Float padding, Bool usezero, Bool useDoublePrec)
118 0 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
119 0 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
120 0 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
121 0 : usePut2_p(false), machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
122 : {
123 0 : mTangent_p=mTangent;
124 0 : tangentSpecified_p=true;
125 0 : useDoubleGrid_p=useDoublePrec;
126 : // peek=NULL;
127 0 : }
128 :
129 6 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
130 : MPosition mLocation, MDirection mTangent, Float padding,
131 6 : Bool usezero, Bool useDoublePrec)
132 6 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
133 6 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
134 6 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
135 12 : usePut2_p(false),machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
136 : {
137 6 : mLocation_p=mLocation;
138 6 : mTangent_p=mTangent;
139 6 : tangentSpecified_p=true;
140 6 : useDoubleGrid_p=useDoublePrec;
141 : // peek=NULL;
142 6 : }
143 :
144 0 : GridFT::GridFT(const RecordInterface& stateRec)
145 0 : : FTMachine()
146 : {
147 : // Construct from the input state record
148 0 : String error;
149 0 : if (!fromRecord(error, stateRec)) {
150 0 : throw (AipsError("Failed to create gridder: " + error));
151 : };
152 0 : timemass_p=0.0;
153 0 : timegrid_p=0.0;
154 0 : timedegrid_p=0.0;
155 : // peek=NULL;
156 0 : }
157 :
158 : //----------------------------------------------------------------------
159 267 : GridFT& GridFT::operator=(const GridFT& other)
160 : {
161 267 : if(this!=&other) {
162 : //Do the base parameters
163 267 : FTMachine::operator=(other);
164 :
165 : //private params
166 267 : imageCache=other.imageCache;
167 267 : cachesize=other.cachesize;
168 267 : tilesize=other.tilesize;
169 267 : convType=other.convType;
170 267 : convType.upcase();
171 267 : uvScale.resize();
172 267 : uvOffset.resize();
173 267 : uvScale.assign(other.uvScale);
174 267 : uvOffset.assign(other.uvOffset);
175 267 : if(other.gridder==0)
176 267 : gridder=0;
177 : else{
178 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
179 0 : uvScale, uvOffset,
180 0 : convType);
181 : }
182 267 : convFunc_p.resize();
183 267 : convSupport_p=other.convSupport_p;
184 267 : convSampling_p=other.convSampling_p;
185 267 : convFunc_p=other.convFunc_p;
186 267 : isTiled=other.isTiled;
187 : //lattice=other.lattice;
188 267 : lattice.reset( );
189 267 : tilesize=other.tilesize;
190 267 : arrayLattice.reset( );
191 267 : maxAbsData=other.maxAbsData;
192 267 : centerLoc=other.centerLoc;
193 267 : offsetLoc=other.offsetLoc;
194 267 : padding_p=other.padding_p;
195 267 : usezero_p=other.usezero_p;
196 267 : noPadding_p=other.noPadding_p;
197 267 : machineName_p="GridFT";
198 267 : timemass_p=0.0;
199 267 : timegrid_p=0.0;
200 : // peek = other.peek;
201 : };
202 267 : return *this;
203 : };
204 :
205 : //----------------------------------------------------------------------
206 267 : GridFT::GridFT(const GridFT& other) : FTMachine(), machineName_p("GridFT")
207 : {
208 267 : operator=(other);
209 267 : }
210 :
211 : //-----------------------------------------------------------------------
212 267 : FTMachine* GridFT::cloneFTM(){
213 267 : return new GridFT(*this);
214 : }
215 :
216 : //===============================
217 2191 : Long GridFT::estimateRAM(const CountedPtr<SIImageStore>& imstor){
218 2191 : Long mem=FTMachine::estimateRAM(imstor);
219 2191 : mem=mem*padding_p*padding_p;
220 2191 : return mem;
221 : }
222 :
223 :
224 : //----------------------------------------------------------------------
225 2948 : void GridFT::init() {
226 :
227 2948 : logIO() << LogOrigin("GridFT", "init") << LogIO::NORMAL;
228 :
229 2948 : 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 2948 : isTiled=false;
250 2948 : if(!noPadding_p && padding_p > 1.01){
251 2926 : CompositeNumber cn(uInt(image->shape()(0)*2));
252 2926 : nx = cn.nextLargerEven(Int(padding_p*Float(image->shape()(0))-0.5));
253 2926 : ny = cn.nextLargerEven(Int(padding_p*Float(image->shape()(1))-0.5));
254 2926 : }
255 : else{
256 22 : nx = image->shape()(0);
257 22 : ny = image->shape()(1);
258 : }
259 2948 : npol = image->shape()(2);
260 2948 : nchan = image->shape()(3);
261 : // }
262 :
263 2948 : sumWeight.resize(npol, nchan);
264 :
265 2948 : uvScale.resize(2);
266 2948 : uvScale(0)=(Float(nx)*image->coordinates().increment()(0));
267 2948 : uvScale(1)=(Float(ny)*image->coordinates().increment()(1));
268 2948 : uvOffset.resize(2);
269 2948 : uvOffset(0)=nx/2;
270 2948 : uvOffset(1)=ny/2;
271 :
272 : // Now set up the gridder. The possibilities are BOX and SF
273 2948 : if(gridder) delete gridder; gridder=0;
274 2948 : convType.upcase();
275 5896 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
276 2948 : uvScale, uvOffset,
277 2948 : convType);
278 :
279 2948 : convFunc_p.resize();
280 2948 : convSupport_p=gridder->cSupport()(0);
281 2948 : convSampling_p=gridder->cSampling();
282 2948 : 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 2948 : if(imageCache) delete imageCache; imageCache=0;
288 :
289 2948 : 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 2948 : }
310 :
311 : // This is nasty, we should use CountedPointers here.
312 5650 : GridFT::~GridFT() {
313 2825 : if(imageCache) delete imageCache; imageCache=0;
314 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
315 2825 : if(gridder) delete gridder; gridder=0;
316 5650 : }
317 :
318 : // Initialize for a transform from the Sky domain. This means that
319 : // we grid-correct, and FFT the image
320 :
321 785 : void GridFT::initializeToVis(ImageInterface<Complex>& iimage,
322 : const vi::VisBuffer2& vb)
323 : {
324 785 : image=&iimage;
325 785 : toVis_p=true;
326 :
327 785 : ok();
328 :
329 785 : init();
330 : // peek->reset();
331 : // Initialize the maps for polarization and channel. These maps
332 : // translate visibility indices into image indices
333 785 : initMaps(vb);
334 : // Need to reset nx, ny for padding
335 : // Padding is possible only for non-tiled processing
336 :
337 :
338 : //cerr << "initialize to vis nx" << nx << " " << ny << " " << npol << " " << nchan << endl;
339 :
340 : //cerr << "image shape " << image->shape() << endl;
341 :
342 : // If we are memory-based then read the image in and create an
343 : // ArrayLattice otherwise just use the PagedImage
344 : /*if(isTiled) {
345 : lattice=std::shared_ptr<Lattice<Complex> >(image, false);
346 : }
347 : else {
348 :
349 : }*/
350 785 : prepGridForDegrid();
351 :
352 : //AlwaysAssert(lattice, AipsError);
353 :
354 :
355 :
356 785 : }
357 :
358 785 : void GridFT::prepGridForDegrid(){
359 785 : IPosition gridShape(4, nx, ny, npol, nchan);
360 785 : griddedData.resize(gridShape);
361 : //griddedData can be a reference of image data...if not using model col
362 : //hence using an undocumented feature of resize that if
363 : //the size is the same as old data it is not changed.
364 : //if(!usePut2_p) griddedData.set(0);
365 785 : griddedData.set(Complex(0.0));
366 :
367 785 : IPosition stride(4, 1);
368 1570 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
369 1570 : IPosition trc(blc+image->shape()-stride);
370 :
371 785 : IPosition start(4, 0);
372 785 : griddedData(blc, trc) = image->getSlice(start, image->shape());
373 :
374 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
375 785 : arrayLattice.reset( new ArrayLattice<Complex>(griddedData) );
376 785 : lattice=arrayLattice;
377 : //logIO() << LogIO::DEBUGGING
378 : // << "Starting grid correction and FFT of image" << LogIO::POST;
379 :
380 : // Do the Grid-correction.
381 : {
382 785 : Vector<Complex> correction(nx);
383 785 : correction=Complex(1.0, 0.0);
384 : // Do the Grid-correction
385 785 : IPosition cursorShape(4, nx, 1, 1, 1);
386 785 : IPosition axisPath(4, 0, 1, 2, 3);
387 785 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
388 785 : LatticeIterator<Complex> lix(*lattice, lsx);
389 555009 : for(lix.reset();!lix.atEnd();lix++) {
390 554224 : gridder->correctX1D(correction, lix.position()(1));
391 554224 : lix.rwVectorCursor()/=correction;
392 : }
393 785 : }
394 785 : image->clearCache();
395 : // Now do the FFT2D in place
396 : //LatticeFFT::cfft2d(*lattice);
397 785 : ft_p.c2cFFT(*lattice);
398 : //logIO() << LogIO::DEBUGGING
399 : // << "Finished grid correction and FFT of image" << LogIO::POST;
400 :
401 785 : }
402 :
403 :
404 785 : void GridFT::finalizeToVis()
405 : {
406 :
407 785 : logIO() << LogOrigin("GridFT", "finalizeToVis") << LogIO::NORMAL;
408 785 : logIO() <<LogIO::NORMAL2<< "Time to degrid " << timedegrid_p <<LogIO::POST;
409 785 : timedegrid_p=0.0;
410 :
411 785 : if(arrayLattice) arrayLattice=nullptr;
412 785 : if(lattice) lattice=nullptr;
413 785 : griddedData.resize();
414 785 : if(isTiled) {
415 :
416 :
417 :
418 0 : AlwaysAssert(imageCache, AipsError);
419 0 : AlwaysAssert(image, AipsError);
420 0 : ostringstream o;
421 0 : imageCache->flush();
422 0 : imageCache->showCacheStatistics(o);
423 0 : logIO() << o.str() << LogIO::POST;
424 0 : }
425 785 : }
426 :
427 :
428 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
429 : // grid.
430 2163 : void GridFT::initializeToSky(ImageInterface<Complex>& iimage,
431 : Matrix<Float>& weight, const vi::VisBuffer2& vb)
432 : {
433 : // image always points to the image
434 2163 : image=&iimage;
435 2163 : toVis_p=false;
436 :
437 2163 : init();
438 :
439 : // Initialize the maps for polarization and channel. These maps
440 : // translate visibility indices into image indices
441 2163 : initMaps(vb);
442 :
443 :
444 :
445 2163 : sumWeight=0.0;
446 2163 : weight.resize(sumWeight.shape());
447 2163 : weight=0.0;
448 :
449 :
450 2163 : IPosition gridShape(4, nx, ny, npol, nchan);
451 2163 : if( !useDoubleGrid_p )
452 : {
453 22 : griddedData.resize(gridShape);
454 22 : griddedData=Complex(0.0);
455 : }
456 : else{
457 :
458 2141 : griddedData2.resize(gridShape);
459 2141 : griddedData2=DComplex(0.0);
460 : }
461 2163 : image->clearCache();
462 : //iimage.get(griddedData, false);
463 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
464 :
465 :
466 : //AlwaysAssert(lattice, AipsError);
467 2163 : }
468 :
469 :
470 :
471 2162 : 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 2162 : logIO() << LogOrigin("GridFT", "finalizeToSky") << LogIO::NORMAL;
477 2162 : logIO()<<LogIO::NORMAL2 <<"Time to massage data " << timemass_p << LogIO::POST;
478 2162 : logIO()<< LogIO::NORMAL2 <<"Time to grid data " << timegrid_p << LogIO::POST;
479 2162 : timemass_p=0.0;
480 2162 : timegrid_p=0.0;
481 2162 : 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 2162 : }
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 524524 : void GridFT::put(const vi::VisBuffer2& vb, Int row, Bool dopsf,
684 : FTMachine::Type type)
685 : {
686 :
687 524524 : 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 524523 : 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 524523 : if(max(chanMap)==-1)
707 : {
708 : // peek->reset();
709 136 : return;
710 : }
711 :
712 524387 : Timer tim;
713 524387 : tim.mark();
714 :
715 : //const Matrix<Float> *imagingweight;
716 : //imagingweight=&(vb.imagingWeight());
717 524387 : Matrix<Float> imagingweight;
718 524387 : getImagingWeight(imagingweight, vb);
719 :
720 524387 : if(dopsf) {type=FTMachine::PSF;}
721 :
722 524387 : Cube<Complex> data;
723 : //Fortran gridder need the flag as ints
724 524387 : Cube<Int> flags;
725 524387 : Matrix<Float> elWeight;
726 524387 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
727 :
728 :
729 : Bool iswgtCopy;
730 : const Float *wgtStorage;
731 524387 : wgtStorage=elWeight.getStorage(iswgtCopy);
732 :
733 : Bool isCopy;
734 : const Complex *datStorage;
735 524387 : if(!dopsf)
736 320866 : datStorage=data.getStorage(isCopy);
737 : else
738 203521 : datStorage=0;
739 : // If row is -1 then we pass through all rows
740 : Int startRow, endRow, nRow;
741 524387 : if (row==-1) {
742 524387 : nRow=vb.nRows();
743 524387 : startRow=0;
744 524387 : 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 524387 : const IPosition &fs=flags.shape();
763 524387 : std::vector<Int> s(fs.begin(),fs.end());
764 524387 : Int nvispol=s[0];
765 524387 : Int nvischan=s[1];
766 524387 : Int nvisrow=s[2];
767 524387 : Matrix<Double> uvw(negateUV(vb));
768 :
769 524387 : Vector<Double> dphase(vb.nRows());
770 524387 : Cube<Int> loc(2, nvischan, nRow);
771 524387 : Cube<Int> off(2, nvischan, nRow);
772 524387 : Matrix<Complex> phasor(nvischan, nRow);
773 524387 : Int csamp=convSampling_p;
774 524387 : dphase=0.0;
775 :
776 524387 : timemass_p +=tim.real();
777 524387 : tim.mark();
778 524387 : rotateUVW(uvw, dphase, vb);
779 524387 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
780 : Bool delphase;
781 524387 : Complex * phasorstor=phasor.getStorage(delphase);
782 524387 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
783 524387 : const Double * scalestor=uvScale.getStorage(del);
784 524387 : const Double * offsetstor=uvOffset.getStorage(del);
785 524387 : const Double* uvstor= uvw.getStorage(del);
786 524387 : Int * locstor=loc.getStorage(del);
787 524387 : Int * offstor=off.getStorage(del);
788 524387 : const Double *dpstor=dphase.getStorage(del);
789 : //Vector<Double> f1=interpVisFreq_p.copy();
790 524387 : Int nvchan=nvischan;
791 : Int irow;
792 524387 : Double cinv=Double(1.0)/C::c;
793 524387 : Int dow=0;
794 524387 : Int nth=1;
795 : #ifdef _OPENMP
796 524387 : if(numthreads_p >0){
797 0 : nth=min(numthreads_p, omp_get_max_threads());
798 : }
799 : else{
800 524387 : nth= omp_get_max_threads();
801 : }
802 : //nth=min(4,nth);
803 : #endif
804 :
805 :
806 524387 : #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 524387 : Int idopsf=0;
819 524387 : 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 524387 : Vector<Int> rowFlags(vb.nRows());
831 524387 : rowFlags=0;
832 524387 : rowFlags(vb.flagRow())=true;
833 524387 : if(!usezero_p) {
834 177502252 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
835 176983879 : 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 524387 : Int csupp=convSupport_p;
843 :
844 524387 : const Double * convfuncstor=(convFunc_p).getStorage(del);
845 :
846 : // cerr <<"Poffset " << min(off) << " " << max(off) << " length " << gridder->cFunction().shape() << endl;
847 :
848 :
849 :
850 524387 : ixsub=1;
851 524387 : iysub=1;
852 524387 : if (nth >4){
853 244106 : ixsub=8;
854 244106 : iysub=8;
855 : }
856 280281 : else if(nth >1) {
857 280281 : ixsub=2;
858 280281 : iysub=2;
859 : }
860 :
861 :
862 524387 : Int rbeg=startRow+1;
863 524387 : Int rend=endRow+1;
864 524387 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
865 17268295 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
866 16743908 : sumwgt[icounter].resize(sumWeight.shape());
867 16743908 : sumwgt[icounter].set(0.0);
868 : }
869 524387 : if(doneThreadPartition_p < 0){
870 1198 : xsect_p.resize(ixsub*iysub);
871 1198 : ysect_p.resize(ixsub*iysub);
872 1198 : nxsect_p.resize(ixsub*iysub);
873 1198 : nysect_p.resize(ixsub*iysub);
874 37130 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
875 35932 : findGridSector(nx, ny, ixsub, iysub, 0, 0, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
876 : }
877 : }
878 524387 : Vector<Int> xsect, ysect, nxsect, nysect;
879 524387 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
880 524387 : const Int* pmapstor=polMap.getStorage(del);
881 524387 : const Int *cmapstor=chanMap.getStorage(del);
882 524387 : Int nc=nchan;
883 524387 : Int np=npol;
884 524387 : Int nxp=nx;
885 524387 : Int nyp=ny;
886 524387 : const Int * flagstor=flags.getStorage(del);
887 524387 : const Int * rowflagstor=rowFlags.getStorage(del);
888 : ////////////////////////
889 :
890 : Bool gridcopy;
891 524387 : if(useDoubleGrid_p){
892 518373 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
893 518373 : #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 17209425 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
929 16691052 : sumWeight=sumWeight+sumwgt[icounter];
930 : }
931 : //phasor.putStorage(phasorstor, delphase);
932 518373 : griddedData2.putStorage(gridstor, gridcopy);
933 518373 : if(dopsf && (nth >4))
934 99298 : tweakGridSector(nx, ny, ixsub, iysub);
935 : }
936 : else{
937 6014 : Complex *gridstor=griddedData.getStorage(gridcopy);
938 6014 : #pragma omp parallel default(none) private(icounter, del) firstprivate(idopsf, datStorage, wgtStorage, flagstor, rowflagstor, convfuncstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc,ixsub, iysub, rend, rbeg, csamp, csupp, nvispol, nvischan, nvisrow, phasorstor, locstor, offstor, xsect, ysect, nxsect, nysect) shared(sumwgt) num_threads(ixsub*iysub)
939 : {
940 : //cerr << "numthreads " << omp_get_num_threads() << endl;
941 : #pragma omp for schedule(dynamic)
942 :
943 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
944 : //cerr << "thread id " << omp_get_thread_num() << endl;
945 : Int x0=xsect(icounter);
946 : Int y0=ysect(icounter);
947 : Int nxsub=nxsect(icounter);
948 : Int nysub=nysect(icounter);
949 :
950 :
951 :
952 : //cerr << "x0 " << x0 << " y0 " << y0 << " nxsub " << nxsub << " nysub " << nysub << endl;
953 : sectggrids(datStorage,
954 : &nvispol,
955 : &nvischan,
956 : &idopsf,
957 : flagstor,
958 : rowflagstor,
959 : wgtStorage,
960 : &nvisrow,
961 : gridstor,
962 : &nxp,
963 : &nyp,
964 : &np,
965 : &nc,
966 : &csupp,
967 : &csamp,
968 : convfuncstor,
969 : cmapstor,
970 : pmapstor,
971 : (sumwgt[icounter]).getStorage(del),
972 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
973 : phasorstor);
974 : }//end for
975 : }// end pragma parallel
976 58870 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
977 52856 : sumWeight=sumWeight+sumwgt[icounter];
978 : }
979 :
980 6014 : griddedData.putStorage(gridstor, gridcopy);
981 6014 : if(dopsf && (nth > 4))
982 480 : tweakGridSector(nx, ny, ixsub, iysub);
983 : }
984 : // cerr << "sunweight " << sumWeight << endl;
985 :
986 524387 : timegrid_p+=tim.real();
987 :
988 524387 : if(!dopsf)
989 320866 : data.freeStorage(datStorage, isCopy);
990 524387 : elWeight.freeStorage(wgtStorage,iswgtCopy);
991 :
992 : // peek->reset();
993 524387 : }
994 :
995 22 : void GridFT::modifyConvFunc(const Vector<Double>& convFunc, Int convSupport, Int convSampling){
996 22 : convFunc_p.resize();
997 22 : convFunc_p=convFunc;
998 22 : convSupport_p=convSupport;
999 22 : convSampling_p=convSampling;
1000 :
1001 22 : }
1002 :
1003 164612 : void GridFT::get(vi::VisBuffer2& vb, Int row)
1004 : {
1005 :
1006 164612 : gridOk(convSupport_p);
1007 : // If row is -1 then we pass through all rows
1008 : Int startRow, endRow, nRow;
1009 164612 : if (row < 0) {
1010 164612 : nRow=vb.nRows();
1011 164612 : startRow=0;
1012 164612 : endRow=nRow-1;
1013 : //unnecessary zeroing
1014 : //vb.modelVisCube()=Complex(0.0,0.0);
1015 : } else {
1016 0 : nRow=1;
1017 0 : startRow=row;
1018 0 : endRow=row;
1019 : //unnecessary zeroing
1020 : //vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
1021 : }
1022 :
1023 :
1024 : ///Channel matching for the actual spectral window of buffer
1025 164612 : matchChannel(vb);
1026 :
1027 :
1028 : //cerr << "GETchanMap " << chanMap << endl;
1029 : //No point in reading data if its not matching in frequency
1030 164612 : if(max(chanMap)==-1)
1031 58 : return;
1032 :
1033 :
1034 :
1035 : // Get the uvws in a form that Fortran can use
1036 164554 : Matrix<Double> uvw(negateUV(vb));
1037 164554 : Vector<Double> dphase(vb.nRows());
1038 164554 : dphase=0.0;
1039 164554 : rotateUVW(uvw, dphase, vb);
1040 164554 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1041 :
1042 :
1043 :
1044 : //Check if ms has changed then cache new spw and chan selection
1045 : //if(vb.isNewMs())
1046 : // matchAllSpwChans(vb);
1047 :
1048 :
1049 : //Here we redo the match or use previous match
1050 :
1051 :
1052 164554 : Cube<Complex> data;
1053 164554 : Cube<Int> flags;
1054 164554 : getInterpolateArrays(vb, data, flags);
1055 :
1056 : //cerr << "max min data " << max(griddedData) << " " << min(griddedData) << " flags " << max(flags) << " " << min(flags) << endl;
1057 :
1058 : // IPosition s(data.shape());
1059 164554 : Int nvp=data.shape()(0);
1060 164554 : Int nvc=data.shape()(1);
1061 164554 : Int nvisrow=data.shape()(2);
1062 :
1063 :
1064 : //cerr << "get flags " << min(flags) << " " << max(flags) << endl;
1065 : Complex *datStorage;
1066 : Bool isCopy;
1067 164554 : datStorage=data.getStorage(isCopy);
1068 :
1069 : ///
1070 164554 : Cube<Int> loc(2, nvc, nvisrow);
1071 164554 : Cube<Int> off(2, nvc, nRow);
1072 164554 : Matrix<Complex> phasor(nvc, nRow);
1073 164554 : Int csamp=convSampling_p;
1074 : Bool delphase;
1075 : Bool del;
1076 164554 : Complex * phasorstor=phasor.getStorage(delphase);
1077 164554 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1078 164554 : const Double * scalestor=uvScale.getStorage(del);
1079 164554 : const Double * offsetstor=uvOffset.getStorage(del);
1080 164554 : const Double* uvstor= uvw.getStorage(del);
1081 164554 : Int * locstor=loc.getStorage(del);
1082 164554 : Int * offstor=off.getStorage(del);
1083 164554 : const Double *dpstor=dphase.getStorage(del);
1084 : //Vector<Double> f1=interpVisFreq_p.copy();
1085 164554 : Int nvchan=nvc;
1086 : Int irow;
1087 164554 : Double cinv=Double(1.0)/C::c;
1088 164554 : Int dow=0;
1089 164554 : Int nth=1;
1090 : #ifdef _OPENMP
1091 164554 : if(numthreads_p >0){
1092 0 : nth=min(numthreads_p, omp_get_max_threads());
1093 : }
1094 : else{
1095 164554 : nth= omp_get_max_threads();
1096 : }
1097 : //nth=min(4,nth);
1098 : #endif
1099 :
1100 :
1101 164554 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvchan, scalestor, offsetstor, csamp, phasorstor, uvstor, locstor, offstor, dpstor, cinv, dow) shared(startRow, endRow) num_threads(nth)
1102 :
1103 : {
1104 : #pragma omp for schedule(dynamic)
1105 : for (irow=startRow; irow<=endRow; ++irow){
1106 : //locateuvw(uvstor,dpstor, visfreqstor, nvchan, scalestor, offsetstor, csamp,
1107 : // locstor,
1108 : // offstor, phasorstor, irow);
1109 :
1110 : locuvw(uvstor, dpstor, visfreqstor, &nvchan, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
1111 : }
1112 :
1113 : }//end pragma parallel
1114 :
1115 164554 : Int rbeg=startRow+1;
1116 164554 : Int rend=endRow+1;
1117 :
1118 :
1119 164554 : Vector<Int> rowFlags(vb.nRows());
1120 164554 : rowFlags=0;
1121 164554 : rowFlags(vb.flagRow())=true;
1122 : //cerr << "rowFlags " << rowFlags << endl;
1123 164554 : if(!usezero_p) {
1124 56762618 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1125 56598064 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
1126 : }
1127 : }
1128 :
1129 :
1130 : //cerr <<"offset " << min(off) << " " <<max(off) << " length " << gridder->cFunction().shape() << endl;
1131 :
1132 : {
1133 : Bool delgrid;
1134 164554 : const Complex* gridstor=griddedData.getStorage(delgrid);
1135 164554 : const Double * convfuncstor=(convFunc_p).getStorage(del);
1136 :
1137 164554 : const Int* pmapstor=polMap.getStorage(del);
1138 164554 : const Int *cmapstor=chanMap.getStorage(del);
1139 164554 : Int nc=nchan;
1140 164554 : Int np=npol;
1141 164554 : Int nxp=nx;
1142 164554 : Int nyp=ny;
1143 164554 : Int csupp=convSupport_p;
1144 164554 : const Int * flagstor=flags.getStorage(del);
1145 164554 : const Int * rowflagstor=rowFlags.getStorage(del);
1146 :
1147 :
1148 164554 : Int npart=nth;
1149 :
1150 :
1151 164554 : Int ix=0;
1152 164554 : #pragma omp parallel default(none) private(ix, rbeg, rend) firstprivate(datStorage, flagstor, rowflagstor, convfuncstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc, csamp, csupp, nvp, nvc, nvisrow, phasorstor, locstor, offstor) shared(npart) num_threads(npart)
1153 : {
1154 : #pragma omp for schedule(dynamic)
1155 : for (ix=0; ix< npart; ++ix){
1156 : rbeg=ix*(nvisrow/npart)+1;
1157 : rend=(ix != (npart-1)) ? (rbeg+(nvisrow/npart)-1) : (rbeg+(nvisrow/npart)+nvisrow%npart-1) ;
1158 : //cerr << "rbeg " << rbeg << " rend " << rend << " " << nvisrow << endl;
1159 :
1160 : sectdgrid(datStorage,
1161 : &nvp,
1162 : &nvc,
1163 : flagstor,
1164 : rowflagstor,
1165 : &nvisrow,
1166 : gridstor,
1167 : &nxp,
1168 : &nyp,
1169 : &np,
1170 : &nc,
1171 : &csupp,
1172 : &csamp,
1173 : convfuncstor,
1174 : cmapstor,
1175 : pmapstor,
1176 : &rbeg, &rend, locstor, offstor, phasorstor);
1177 : }//end pragma parallel
1178 : }
1179 164554 : data.putStorage(datStorage, isCopy);
1180 164554 : griddedData.freeStorage(gridstor, delgrid);
1181 : //cerr << "Get min max " << min(data) << " " << max(data) << endl;
1182 : }
1183 164554 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1184 :
1185 164554 : }
1186 :
1187 :
1188 :
1189 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1190 : // return the resulting image
1191 2140 : ImageInterface<Complex>& GridFT::getImage(Matrix<Float>& weights, Bool normalize)
1192 : {
1193 : //AlwaysAssert(lattice, AipsError);
1194 2140 : AlwaysAssert(gridder, AipsError);
1195 2140 : AlwaysAssert(image, AipsError);
1196 2140 : logIO() << LogOrigin("GridFT", "getImage") << LogIO::NORMAL;
1197 :
1198 2140 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
1199 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
1200 : }
1201 2140 : weights.resize(sumWeight.shape());
1202 :
1203 2140 : convertArray(weights, sumWeight);
1204 : // If the weights are all zero then we cannot normalize
1205 : // otherwise we don't care.
1206 2140 : 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 2140 : if(useDoubleGrid_p)
1233 : {
1234 2140 : ArrayLattice<DComplex> darrayLattice(griddedData2);
1235 : //LatticeFFT::cfft2d(darrayLattice,false);
1236 2140 : ft_p.c2cFFT(darrayLattice, False);
1237 2140 : griddedData.resize(griddedData2.shape());
1238 2140 : convertArray(griddedData, griddedData2);
1239 :
1240 2140 : SynthesisUtilMethods::getResource("mem peak in getImage");
1241 :
1242 : //Don't need the double-prec grid anymore...
1243 2140 : griddedData2.resize();
1244 2140 : arrayLattice.reset( new ArrayLattice<Complex>(griddedData) );
1245 2140 : lattice=arrayLattice;
1246 2140 : }
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 2140 : Int inx = lattice->shape()(0);
1259 2140 : Int iny = lattice->shape()(1);
1260 2140 : Vector<Complex> correction(inx);
1261 2140 : correction=Complex(1.0, 0.0);
1262 : // Do the Grid-correction
1263 2140 : IPosition cursorShape(4, inx, 1, 1, 1);
1264 2140 : IPosition axisPath(4, 0, 1, 2, 3);
1265 2140 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1266 2140 : LatticeIterator<Complex> lix(*lattice, lsx);
1267 1582588 : for(lix.reset();!lix.atEnd();lix++) {
1268 1580448 : Int pol=lix.position()(2);
1269 1580448 : Int chan=lix.position()(3);
1270 1580448 : if(weights(pol, chan)!=0.0) {
1271 1523400 : gridder->correctX1D(correction, lix.position()(1));
1272 1523400 : lix.rwVectorCursor()/=correction;
1273 1523400 : if(normalize) {
1274 0 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
1275 0 : lix.rwCursor()*=rnorm;
1276 : }
1277 : else {
1278 1523400 : Complex rnorm(Float(inx)*Float(iny));
1279 1523400 : lix.rwCursor()*=rnorm;
1280 : }
1281 : }
1282 : else {
1283 57048 : lix.woCursor()=0.0;
1284 : }
1285 : }
1286 2140 : }
1287 2140 : LatticeLocker lock1 (*(image), FileLocker::Write);
1288 2140 : if(!isTiled) {
1289 : // Check the section from the image BEFORE converting to a lattice
1290 4280 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1291 2140 : IPosition stride(4, 1);
1292 4280 : IPosition trc(blc+image->shape()-stride);
1293 : // Do the copy
1294 2140 : IPosition start(4, 0);
1295 :
1296 2140 : image->put(griddedData(blc, trc));
1297 2140 : }
1298 2140 : }
1299 2140 : image->flush();
1300 2140 : image->clearCache();
1301 :
1302 2140 : if(arrayLattice) arrayLattice=nullptr;
1303 2140 : if(lattice) lattice=nullptr;
1304 2140 : griddedData.resize();
1305 2140 : return *image;
1306 : }
1307 :
1308 : // Get weight image
1309 0 : void GridFT::getWeightImage(ImageInterface<Float>& weightImage, Matrix<Float>& weights)
1310 : {
1311 :
1312 0 : logIO() << LogOrigin("GridFT", "getWeightImage") << LogIO::NORMAL;
1313 :
1314 0 : weights.resize(sumWeight.shape());
1315 0 : convertArray(weights,sumWeight);
1316 :
1317 0 : const IPosition latticeShape = weightImage.shape();
1318 :
1319 0 : Int nx=latticeShape(0);
1320 0 : Int ny=latticeShape(1);
1321 :
1322 0 : IPosition loc(2, 0);
1323 0 : IPosition cursorShape(4, nx, ny, 1, 1);
1324 0 : IPosition axisPath(4, 0, 1, 2, 3);
1325 0 : LatticeStepper lsx(latticeShape, cursorShape, axisPath);
1326 0 : LatticeIterator<Float> lix(weightImage, lsx);
1327 0 : for(lix.reset();!lix.atEnd();lix++) {
1328 0 : Int pol=lix.position()(2);
1329 0 : Int chan=lix.position()(3);
1330 0 : lix.rwCursor()=weights(pol,chan);
1331 : }
1332 0 : }
1333 :
1334 45 : Bool GridFT::toRecord(String& error,
1335 : RecordInterface& outRec, Bool withImage, const String diskimage)
1336 : {
1337 :
1338 :
1339 :
1340 : // Save the current GridFT object to an output state record
1341 45 : Bool retval = true;
1342 45 : Float elpadd=padding_p;
1343 45 : if(toVis_p && withImage)
1344 19 : elpadd=1.0;
1345 : //save the base class variables
1346 45 : if(!FTMachine::toRecord(error, outRec, withImage, diskimage))
1347 0 : return false;
1348 :
1349 : //a call to init will redo imagecache and gridder
1350 : // so no need to save these
1351 :
1352 45 : outRec.define("cachesize", Int64(cachesize));
1353 45 : outRec.define("tilesize", tilesize);
1354 45 : outRec.define("convtype", convType);
1355 45 : outRec.define("uvscale", uvScale);
1356 45 : outRec.define("uvoffset", uvOffset);
1357 45 : outRec.define("istiled", isTiled);
1358 45 : outRec.define("maxabsdata", maxAbsData);
1359 45 : Vector<Int> center_loc(4), offset_loc(4);
1360 225 : for (Int k=0; k<4 ; k++){
1361 180 : center_loc(k)=centerLoc(k);
1362 180 : offset_loc(k)=offsetLoc(k);
1363 : }
1364 45 : outRec.define("centerloc", center_loc);
1365 45 : outRec.define("offsetloc", offset_loc);
1366 45 : outRec.define("padding", elpadd);
1367 45 : outRec.define("usezero", usezero_p);
1368 45 : outRec.define("nopadding", noPadding_p);
1369 45 : return retval;
1370 45 : }
1371 :
1372 0 : Bool GridFT::fromRecord(String& error,
1373 : const RecordInterface& inRec)
1374 : {
1375 0 : Bool retval = true;
1376 0 : if(!FTMachine::fromRecord(error, inRec))
1377 0 : return false;
1378 0 : gridder=0; imageCache=0; lattice.reset( ); arrayLattice.reset( );
1379 : //For some reason int64 seems to be getting lost...
1380 : //this is not an important parameter...so filing a jira
1381 0 : if(inRec.isDefined("cachesize"))
1382 0 : cachesize=inRec.asInt64("cachesize");
1383 : else
1384 0 : cachesize=1000000;
1385 0 : inRec.get("tilesize", tilesize);
1386 0 : inRec.get("convtype", convType);
1387 0 : inRec.get("uvscale", uvScale);
1388 0 : inRec.get("uvoffset", uvOffset);
1389 0 : inRec.get("istiled", isTiled);
1390 0 : inRec.get("maxabsdata", maxAbsData);
1391 0 : Vector<Int> center_loc(4), offset_loc(4);
1392 0 : inRec.get("centerloc", center_loc);
1393 0 : inRec.get("offsetloc", offset_loc);
1394 0 : uInt ndim4 = 4;
1395 0 : centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2),
1396 0 : center_loc(3));
1397 0 : offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2),
1398 0 : offset_loc(3));
1399 0 : inRec.get("padding", padding_p);
1400 0 : inRec.get("usezero", usezero_p);
1401 0 : inRec.get("nopadding", noPadding_p);
1402 :
1403 0 : machineName_p="GridFT";
1404 : ///setup some of the parameters if there is an image
1405 0 : if(inRec.isDefined("image") || inRec.isDefined("diskimage"))
1406 0 : init();
1407 : /*if(!cmplxImage_p.null()){
1408 : //FTMachine::fromRecord would have recovered the image
1409 : // Might be changing the shape of sumWeight
1410 :
1411 : if(isTiled) {
1412 : lattice=std::shared_ptr<Lattice<Complex> >(image, false);
1413 : }
1414 : else {
1415 : // Make the grid the correct shape and turn it into an array lattice
1416 : // Check the section from the image BEFORE converting to a lattice
1417 : if(!toVis_p){
1418 : IPosition gridShape(4, nx, ny, npol, nchan);
1419 : griddedData.resize(gridShape);
1420 : griddedData=Complex(0.0);
1421 : }
1422 : else{
1423 : prepGridForDegrid();
1424 : }
1425 : }
1426 : };*/
1427 0 : return retval;
1428 0 : }
1429 :
1430 1381615 : void GridFT::ok() {
1431 1381615 : AlwaysAssert(image, AipsError);
1432 1381615 : }
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 4140 : String GridFT::name() const {
1502 :
1503 4140 : return machineName_p;
1504 :
1505 :
1506 : }
1507 :
1508 : }//End of namespace refim
1509 : } //# NAMESPACE CASA - END
1510 :
|
