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