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