Line data Source code
1 : //# GridFT.cc: Implementation of GridFT class
2 : //# Copyright (C) 1997-2014
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 Library 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 Library 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/VisibilityIterator.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 <synthesis/TransformMachines/GridFT.h>
41 : #include <casacore/scimath/Mathematics/RigidVector.h>
42 : #include <msvis/MSVis/StokesVector.h>
43 : #include <synthesis/TransformMachines/StokesImageUtil.h>
44 : #include <synthesis/TransformMachines/Utils.h>
45 : #include <msvis/MSVis/VisBuffer.h>
46 : #include <msvis/MSVis/VisSet.h>
47 : #include <casacore/images/Images/ImageInterface.h>
48 : #include <casacore/images/Images/PagedImage.h>
49 : #include <casacore/casa/Containers/Block.h>
50 : #include <casacore/casa/Containers/Record.h>
51 : #include <casacore/casa/Arrays/ArrayLogical.h>
52 : #include <casacore/casa/Arrays/ArrayMath.h>
53 : #include <casacore/casa/Arrays/Array.h>
54 : #include <casacore/casa/Arrays/MaskedArray.h>
55 : #include <casacore/casa/Arrays/Vector.h>
56 : #include <casacore/casa/Arrays/Matrix.h>
57 : #include <casacore/casa/Arrays/Cube.h>
58 : #include <casacore/casa/Arrays/MatrixIter.h>
59 : #include <casacore/casa/BasicSL/String.h>
60 : #include <casacore/casa/Utilities/Assert.h>
61 : #include <casacore/casa/Exceptions/Error.h>
62 : #include <casacore/lattices/Lattices/ArrayLattice.h>
63 : #include <casacore/measures/Measures/UVWMachine.h>
64 : #include <casacore/lattices/Lattices/SubLattice.h>
65 : #include <casacore/lattices/LRegions/LCBox.h>
66 : #include <casacore/lattices/Lattices/LatticeCache.h>
67 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
68 : #include <casacore/lattices/Lattices/LatticeIterator.h>
69 : #include <casacore/lattices/Lattices/LatticeStepper.h>
70 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
71 : #include <casacore/casa/Utilities/CompositeNumber.h>
72 : #include <casacore/casa/OS/Timer.h>
73 : #include <sstream>
74 : #ifdef _OPENMP
75 : #include <omp.h>
76 : #endif
77 :
78 : using namespace casacore;
79 : namespace casa { //# NAMESPACE CASA - BEGIN
80 : using namespace casacore;
81 : // using namespace casa::async;
82 0 : GridFT::GridFT() : FTMachine(), padding_p(1.0), imageCache(0), cachesize(1000000), tilesize(1000), gridder(0), isTiled(false), convType("SF"),
83 0 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
84 0 : usezero_p(false), noPadding_p(false), usePut2_p(false),
85 0 : machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0) {
86 :
87 0 : }
88 0 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType, Float padding,
89 0 : Bool usezero, Bool useDoublePrec)
90 0 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize), tilesize(itilesize),
91 0 : gridder(0), isTiled(false), convType(iconvType),
92 0 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
93 0 : usezero_p(usezero), noPadding_p(false), usePut2_p(false),
94 0 : machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
95 : {
96 0 : useDoubleGrid_p=useDoublePrec;
97 : // peek=NULL;
98 0 : }
99 :
100 29 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
101 29 : MPosition mLocation, Float padding, Bool usezero, Bool useDoublePrec)
102 29 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
103 29 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
104 29 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
105 58 : usePut2_p(false), machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
106 : {
107 29 : mLocation_p=mLocation;
108 29 : tangentSpecified_p=false;
109 29 : useDoubleGrid_p=useDoublePrec;
110 : // peek=NULL;
111 29 : }
112 :
113 0 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
114 0 : MDirection mTangent, Float padding, Bool usezero, Bool useDoublePrec)
115 0 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
116 0 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
117 0 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
118 0 : usePut2_p(false), machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
119 : {
120 0 : mTangent_p=mTangent;
121 0 : tangentSpecified_p=true;
122 0 : useDoubleGrid_p=useDoublePrec;
123 : // peek=NULL;
124 0 : }
125 :
126 0 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
127 : MPosition mLocation, MDirection mTangent, Float padding,
128 0 : Bool usezero, Bool useDoublePrec)
129 0 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
130 0 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
131 0 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
132 0 : usePut2_p(false),machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
133 : {
134 0 : mLocation_p=mLocation;
135 0 : mTangent_p=mTangent;
136 0 : tangentSpecified_p=true;
137 0 : useDoubleGrid_p=useDoublePrec;
138 : // peek=NULL;
139 0 : }
140 :
141 0 : GridFT::GridFT(const RecordInterface& stateRec)
142 0 : : FTMachine()
143 : {
144 : // Construct from the input state record
145 0 : String error;
146 0 : if (!fromRecord(error, stateRec)) {
147 0 : throw (AipsError("Failed to create gridder: " + error));
148 : };
149 0 : timemass_p=0.0;
150 0 : timegrid_p=0.0;
151 0 : timedegrid_p=0.0;
152 : // peek=NULL;
153 0 : }
154 :
155 : //----------------------------------------------------------------------
156 28 : GridFT& GridFT::operator=(const GridFT& other)
157 : {
158 28 : if(this!=&other) {
159 : //Do the base parameters
160 28 : FTMachine::operator=(other);
161 :
162 : //private params
163 28 : imageCache=other.imageCache;
164 28 : cachesize=other.cachesize;
165 28 : tilesize=other.tilesize;
166 28 : convType=other.convType;
167 28 : convType.upcase();
168 28 : uvScale.resize();
169 28 : uvOffset.resize();
170 28 : uvScale.assign(other.uvScale);
171 28 : uvOffset.assign(other.uvOffset);
172 28 : if(other.gridder==0)
173 28 : gridder=0;
174 : else{
175 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
176 0 : uvScale, uvOffset,
177 0 : convType);
178 : }
179 28 : convFunc_p.resize();
180 28 : convSupport_p=other.convSupport_p;
181 28 : convSampling_p=other.convSampling_p;
182 28 : convFunc_p=other.convFunc_p;
183 28 : isTiled=other.isTiled;
184 : //lattice=other.lattice;
185 28 : lattice=0;
186 28 : tilesize=other.tilesize;
187 28 : arrayLattice=0;
188 28 : maxAbsData=other.maxAbsData;
189 28 : centerLoc=other.centerLoc;
190 28 : offsetLoc=other.offsetLoc;
191 28 : padding_p=other.padding_p;
192 28 : usezero_p=other.usezero_p;
193 28 : noPadding_p=other.noPadding_p;
194 28 : machineName_p="GridFT";
195 28 : timemass_p=0.0;
196 28 : timegrid_p=0.0;
197 : // peek = other.peek;
198 : };
199 28 : return *this;
200 : };
201 :
202 : //----------------------------------------------------------------------
203 28 : GridFT::GridFT(const GridFT& other) : FTMachine(), machineName_p("GridFT")
204 : {
205 28 : operator=(other);
206 28 : }
207 :
208 : //-----------------------------------------------------------------------
209 0 : FTMachine* GridFT::cloneFTM(){
210 0 : return new GridFT(*this);
211 : }
212 :
213 : //----------------------------------------------------------------------
214 55 : void GridFT::init() {
215 :
216 55 : logIO() << LogOrigin("GridFT", "init") << LogIO::NORMAL;
217 :
218 55 : ok();
219 : // if (peek == NULL)
220 : // {
221 : // // peek = new SynthesisAsyncPeek();
222 : // // peek->reset();
223 : // // peek->startThread();
224 : // }
225 :
226 : /* hardwiring isTiled is false
227 : // Padding is possible only for non-tiled processing
228 : if((padding_p*padding_p*image->shape().product())>cachesize) {
229 : isTiled=true;
230 : nx = image->shape()(0);
231 : ny = image->shape()(1);
232 : npol = image->shape()(2);
233 : nchan = image->shape()(3);
234 : }
235 : else {
236 : */
237 : // We are padding.
238 55 : isTiled=false;
239 55 : if(!noPadding_p){
240 55 : CompositeNumber cn(uInt(image->shape()(0)*2));
241 55 : nx = cn.nextLargerEven(Int(padding_p*Float(image->shape()(0))-0.5));
242 55 : ny = cn.nextLargerEven(Int(padding_p*Float(image->shape()(1))-0.5));
243 55 : }
244 : else{
245 0 : nx = image->shape()(0);
246 0 : ny = image->shape()(1);
247 : }
248 55 : npol = image->shape()(2);
249 55 : nchan = image->shape()(3);
250 : // }
251 :
252 55 : sumWeight.resize(npol, nchan);
253 :
254 55 : uvScale.resize(2);
255 55 : uvScale(0)=(Float(nx)*image->coordinates().increment()(0));
256 55 : uvScale(1)=(Float(ny)*image->coordinates().increment()(1));
257 55 : uvOffset.resize(2);
258 55 : uvOffset(0)=nx/2;
259 55 : uvOffset(1)=ny/2;
260 :
261 : // Now set up the gridder. The possibilities are BOX and SF
262 55 : if(gridder) delete gridder; gridder=0;
263 55 : convType.upcase();
264 110 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
265 55 : uvScale, uvOffset,
266 55 : convType);
267 :
268 :
269 55 : convFunc_p.resize();
270 55 : convSupport_p=gridder->cSupport()(0);
271 55 : convSampling_p=gridder->cSampling();
272 55 : convFunc_p=gridder->cFunction();
273 : // Set up image cache needed for gridding. For BOX-car convolution
274 : // we can use non-overlapped tiles. Otherwise we need to use
275 : // overlapped tiles and additive gridding so that only increments
276 : // to a tile are written.
277 55 : if(imageCache) delete imageCache; imageCache=0;
278 :
279 55 : if(isTiled) {
280 0 : Float tileOverlap=0.5;
281 0 : if(convType=="BOX") {
282 0 : tileOverlap=0.0;
283 : }
284 : else {
285 0 : tileOverlap=0.5;
286 0 : tilesize=max(12,tilesize);
287 : }
288 0 : IPosition tileShape=IPosition(4,tilesize,tilesize,npol,nchan);
289 0 : Vector<Float> tileOverlapVec(4);
290 0 : tileOverlapVec=0.0;
291 0 : tileOverlapVec(0)=tileOverlap;
292 0 : tileOverlapVec(1)=tileOverlap;
293 0 : Int tmpCacheVal=static_cast<Int>(cachesize);
294 0 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
295 : tileOverlapVec,
296 0 : (tileOverlap>0.0));
297 :
298 0 : }
299 55 : }
300 :
301 : // This is nasty, we should use CountedPointers here.
302 112 : GridFT::~GridFT() {
303 56 : if(imageCache) delete imageCache; imageCache=0;
304 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
305 56 : if(gridder) delete gridder; gridder=0;
306 112 : }
307 :
308 : // Initialize for a transform from the Sky domain. This means that
309 : // we grid-correct, and FFT the image
310 :
311 40 : void GridFT::initializeToVis(ImageInterface<Complex>& iimage,
312 : const VisBuffer& vb)
313 : {
314 40 : image=&iimage;
315 40 : toVis_p=true;
316 :
317 40 : ok();
318 :
319 40 : init();
320 : // peek->reset();
321 : // Initialize the maps for polarization and channel. These maps
322 : // translate visibility indices into image indices
323 40 : initMaps(vb);
324 :
325 : // Need to reset nx, ny for padding
326 : // Padding is possible only for non-tiled processing
327 :
328 :
329 : //cerr << "initialize to vis nx" << nx << " " << ny << " " << npol << " " << nchan << endl;
330 :
331 : //cerr << "image shape " << image->shape() << endl;
332 :
333 : // If we are memory-based then read the image in and create an
334 : // ArrayLattice otherwise just use the PagedImage
335 : /*if(isTiled) {
336 : lattice=CountedPtr<Lattice<Complex> >(image, false);
337 : }
338 : else {
339 :
340 : }*/
341 40 : prepGridForDegrid();
342 :
343 : //AlwaysAssert(lattice, AipsError);
344 :
345 :
346 :
347 40 : }
348 :
349 40 : void GridFT::prepGridForDegrid(){
350 40 : IPosition gridShape(4, nx, ny, npol, nchan);
351 40 : griddedData.resize(gridShape);
352 : //griddedData can be a reference of image data...if not using model col
353 : //hence using an undocumented feature of resize that if
354 : //the size is the same as old data it is not changed.
355 : //if(!usePut2_p) griddedData.set(0);
356 40 : griddedData.set(Complex(0.0));
357 :
358 40 : IPosition stride(4, 1);
359 80 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
360 80 : IPosition trc(blc+image->shape()-stride);
361 :
362 40 : IPosition start(4, 0);
363 40 : griddedData(blc, trc) = image->getSlice(start, image->shape());
364 :
365 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
366 40 : arrayLattice = new ArrayLattice<Complex>(griddedData);
367 40 : lattice=arrayLattice;
368 : //logIO() << LogIO::DEBUGGING
369 : // << "Starting grid correction and FFT of image" << LogIO::POST;
370 :
371 : // Do the Grid-correction.
372 : {
373 40 : Vector<Complex> correction(nx);
374 40 : correction=Complex(1.0, 0.0);
375 : // Do the Grid-correction
376 40 : IPosition cursorShape(4, nx, 1, 1, 1);
377 40 : IPosition axisPath(4, 0, 1, 2, 3);
378 40 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
379 40 : LatticeIterator<Complex> lix(*lattice, lsx);
380 100040 : for(lix.reset();!lix.atEnd();lix++) {
381 100000 : gridder->correctX1D(correction, lix.position()(1));
382 100000 : lix.rwVectorCursor()/=correction;
383 : }
384 40 : }
385 40 : image->clearCache();
386 : // Now do the FFT2D in place
387 40 : LatticeFFT::cfft2d(*lattice);
388 :
389 : //logIO() << LogIO::DEBUGGING
390 : // << "Finished grid correction and FFT of image" << LogIO::POST;
391 :
392 40 : }
393 :
394 :
395 0 : void GridFT::finalizeToVis()
396 : {
397 :
398 : //cerr <<"Time to degrid " << timedegrid_p << endl;
399 0 : timedegrid_p=0.0;
400 :
401 0 : if(!arrayLattice.null()) arrayLattice=0;
402 0 : if(!lattice.null()) lattice=0;
403 0 : griddedData.resize();
404 :
405 0 : if(isTiled) {
406 :
407 0 : logIO() << LogOrigin("GridFT", "finalizeToVis") << LogIO::NORMAL;
408 :
409 0 : AlwaysAssert(imageCache, AipsError);
410 0 : AlwaysAssert(image, AipsError);
411 0 : ostringstream o;
412 0 : imageCache->flush();
413 0 : imageCache->showCacheStatistics(o);
414 0 : logIO() << o.str() << LogIO::POST;
415 0 : }
416 0 : }
417 :
418 :
419 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
420 : // grid.
421 15 : void GridFT::initializeToSky(ImageInterface<Complex>& iimage,
422 : Matrix<Float>& weight, const VisBuffer& vb)
423 : {
424 : // image always points to the image
425 15 : image=&iimage;
426 15 : toVis_p=false;
427 :
428 15 : init();
429 :
430 : // Initialize the maps for polarization and channel. These maps
431 : // translate visibility indices into image indices
432 15 : initMaps(vb);
433 :
434 :
435 :
436 15 : sumWeight=0.0;
437 15 : weight.resize(sumWeight.shape());
438 15 : weight=0.0;
439 :
440 : // Initialize for in memory or to disk gridding. lattice will
441 : // point to the appropriate Lattice, either the ArrayLattice for
442 : // in memory gridding or to the image for to disk gridding.
443 15 : if(isTiled) {
444 0 : imageCache->flush();
445 0 : image->set(Complex(0.0));
446 0 : lattice=CountedPtr<Lattice<Complex> >(image, false);
447 : }
448 : else {
449 15 : IPosition gridShape(4, nx, ny, npol, nchan);
450 15 : if( !useDoubleGrid_p )
451 : {
452 1 : griddedData.resize(gridShape);
453 1 : griddedData=Complex(0.0);
454 : }
455 15 : if(useDoubleGrid_p){
456 : //griddedData.resize();
457 14 : griddedData2.resize(gridShape);
458 14 : griddedData2=DComplex(0.0);
459 : }
460 15 : image->clearCache();
461 : //iimage.get(griddedData, false);
462 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
463 :
464 15 : }
465 : //AlwaysAssert(lattice, AipsError);
466 15 : }
467 :
468 :
469 :
470 15 : void GridFT::finalizeToSky()
471 : {
472 : //AlwaysAssert(lattice, AipsError);
473 : // Now we flush the cache and report statistics
474 : // For memory based, we don't write anything out yet.
475 : //cerr <<"Time to massage data " << timemass_p << endl;
476 : //cerr <<"Time to grid data " << timegrid_p << endl;
477 15 : timemass_p=0.0;
478 15 : timegrid_p=0.0;
479 15 : if(isTiled) {
480 0 : logIO() << LogOrigin("GridFT", "finalizeToSky") << LogIO::NORMAL;
481 :
482 0 : AlwaysAssert(image, AipsError);
483 0 : AlwaysAssert(imageCache, AipsError);
484 0 : imageCache->flush();
485 0 : ostringstream o;
486 0 : imageCache->showCacheStatistics(o);
487 0 : logIO() << o.str() << LogIO::POST;
488 0 : }
489 : // peek->terminate();
490 : // peek->reset();
491 15 : }
492 :
493 :
494 :
495 0 : Array<Complex>* GridFT::getDataPointer(const IPosition& centerLoc2D,
496 : Bool readonly) {
497 : Array<Complex>* result;
498 : // Is tiled: get tiles and set up offsets
499 0 : centerLoc(0)=centerLoc2D(0);
500 0 : centerLoc(1)=centerLoc2D(1);
501 0 : result=&imageCache->tile(offsetLoc,centerLoc, readonly);
502 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
503 0 : return result;
504 : }
505 :
506 : #define NEED_UNDERSCORES
507 : #if defined(NEED_UNDERSCORES)
508 : #define ggrid ggrid_
509 : #define dgrid dgrid_
510 : #define ggrids ggrids_
511 : #define sectggridd sectggridd_
512 : #define sectggrids sectggrids_
513 : #define sectdgrid sectdgrid_
514 : #define locuvw locuvw_
515 : #endif
516 :
517 : extern "C" {
518 : void locuvw(const Double*, const Double*, const Double*, const Int*, const Double*, const Double*, const Int*,
519 : Int*, Int*, Complex*, const Int*, const Int*, const Double*);
520 : void ggrid(Double*,
521 : Double*,
522 : const Complex*,
523 : Int*,
524 : Int*,
525 : Int*,
526 : const Int*,
527 : const Int*,
528 : const Float*,
529 : Int*,
530 : Int*,
531 : Double*,
532 : Double*,
533 : DComplex*,
534 : Int*,
535 : Int*,
536 : Int *,
537 : Int *,
538 : const Double*,
539 : const Double*,
540 : Int*,
541 : Int*,
542 : Double*,
543 : Int*,
544 : Int*,
545 : Double*);
546 : void sectggridd(const Complex*,
547 : const Int*,
548 : const Int*,
549 : const Int*,
550 : const Int*,
551 : const Int*,
552 : const Float*,
553 : const Int*,
554 : DComplex*,
555 : const Int*,
556 : const Int*,
557 : const Int *,
558 : const Int *,
559 : //support
560 : const Int*,
561 : const Int*,
562 : const Double*,
563 : const Int*,
564 : const Int*,
565 : Double*,
566 : const Int*,
567 : const Int*,
568 : const Int*,
569 : const Int*,
570 : const Int*,
571 : const Int*,
572 : //loc, off, phasor
573 : const Int*,
574 : const Int*,
575 : const Complex*);
576 : //single precision gridder
577 : void sectggrids(const Complex*,
578 : const Int*,
579 : const Int*,
580 : const Int*,
581 : const Int*,
582 : const Int*,
583 : const Float*,
584 : const Int*,
585 : Complex*,
586 : const Int*,
587 : const Int*,
588 : const Int *,
589 : const Int *,
590 : const Int*,
591 : const Int*,
592 : const Double*,
593 : const Int*,
594 : const Int*,
595 : Double*,
596 : const Int*,
597 : const Int*,
598 : const Int*,
599 : const Int*,
600 : const Int*,
601 : const Int*,
602 : //loc, off, phasor
603 : const Int*,
604 : const Int*,
605 : const Complex*);
606 : void ggrids(Double*,
607 : Double*,
608 : const Complex*,
609 : Int*,
610 : Int*,
611 : Int*,
612 : const Int*,
613 : const Int*,
614 : const Float*,
615 : Int*,
616 : Int*,
617 : Double*,
618 : Double*,
619 : Complex*,
620 : Int*,
621 : Int*,
622 : Int *,
623 : Int *,
624 : const Double*,
625 : const Double*,
626 : Int*,
627 : Int*,
628 : Double*,
629 : Int*,
630 : Int*,
631 : Double*);
632 :
633 : void dgrid(Double*,
634 : Double*,
635 : Complex*,
636 : Int*,
637 : Int*,
638 : const Int*,
639 : const Int*,
640 : Int*,
641 : Int*,
642 : Double*,
643 : Double*,
644 : const Complex*,
645 : Int*,
646 : Int*,
647 : Int *,
648 : Int *,
649 : const Double*,
650 : const Double*,
651 : Int*,
652 : Int*,
653 : Double*,
654 : Int*,
655 : Int*);
656 :
657 : void sectdgrid(Complex*,
658 : const Int*,
659 : const Int*,
660 : const Int*,
661 : const Int*,
662 : const Int*,
663 : const Complex*,
664 : const Int*,
665 : const Int*,
666 : const Int *,
667 : const Int *,
668 : //support
669 : const Int*,
670 : const Int*,
671 : const Double*,
672 : const Int*,
673 : const Int*,
674 : //rbeg, rend, loc, off, phasor
675 : const Int*,
676 : const Int*,
677 : const Int*,
678 : const Int*,
679 : const Complex*);
680 : }
681 463 : void GridFT::put(const VisBuffer& vb, Int row, Bool dopsf,
682 : FTMachine::Type type)
683 : {
684 :
685 463 : gridOk(convSupport_p);
686 : // peek->setVBPtr(&vb);
687 :
688 : //Check if ms has changed then cache new spw and chan selection
689 463 : if(vb.newMS())
690 0 : matchAllSpwChans(vb);
691 :
692 : //Here we redo the match or use previous match
693 :
694 : //Channel matching for the actual spectral window of buffer
695 463 : if(doConversion_p[vb.spectralWindow()]){
696 0 : matchChannel(vb.spectralWindow(), vb);
697 : }
698 : else{
699 463 : chanMap.resize();
700 463 : chanMap=multiChanMap_p[vb.spectralWindow()];
701 : }
702 :
703 : //No point in reading data if its not matching in frequency
704 463 : if(max(chanMap)==-1)
705 : {
706 : // peek->reset();
707 0 : return;
708 : }
709 :
710 463 : Timer tim;
711 463 : tim.mark();
712 :
713 : const Matrix<Float> *imagingweight;
714 463 : imagingweight=&(vb.imagingWeight());
715 :
716 463 : if(dopsf) {type=FTMachine::PSF;}
717 :
718 463 : Cube<Complex> data;
719 : //Fortran gridder need the flag as ints
720 463 : Cube<Int> flags;
721 463 : Matrix<Float> elWeight;
722 463 : interpolateFrequencyTogrid(vb, *imagingweight,data, flags, elWeight, type);
723 :
724 :
725 : Bool iswgtCopy;
726 : const Float *wgtStorage;
727 463 : wgtStorage=elWeight.getStorage(iswgtCopy);
728 :
729 : Bool isCopy;
730 : const Complex *datStorage;
731 463 : if(!dopsf)
732 200 : datStorage=data.getStorage(isCopy);
733 : else
734 263 : datStorage=0;
735 : // If row is -1 then we pass through all rows
736 : Int startRow, endRow, nRow;
737 463 : if (row==-1) {
738 463 : nRow=vb.nRow();
739 463 : startRow=0;
740 463 : endRow=nRow-1;
741 : } else {
742 0 : nRow=1;
743 0 : startRow=row;
744 0 : endRow=row;
745 : }
746 :
747 : //cerr << "nRow " << nRow << endl;
748 :
749 :
750 : // Get the uvws in a form that Fortran can use and do that
751 : // necessary phase rotation. On a Pentium Pro 200 MHz
752 : // when null, this step takes about 50us per uvw point. This
753 : // is just barely noticeable for Stokes I continuum and
754 : // irrelevant for other cases.
755 :
756 : Bool del;
757 : // IPosition s(flags.shape());
758 463 : const IPosition &fs=flags.shape();
759 463 : std::vector<Int> s(fs.begin(),fs.end());
760 463 : Int nvispol=s[0];
761 463 : Int nvischan=s[1];
762 463 : Int nvisrow=s[2];
763 463 : Matrix<Double> uvw(3, vb.uvw().nelements());
764 463 : uvw=0.0;
765 463 : Vector<Double> dphase(vb.uvw().nelements());
766 463 : Cube<Int> loc(2, nvischan, nRow);
767 463 : Cube<Int> off(2, nvischan, nRow);
768 463 : Matrix<Complex> phasor(nvischan, nRow);
769 463 : Int csamp=convSampling_p;
770 463 : dphase=0.0;
771 : //NEGATING to correct for an image inversion problem
772 320178 : for (Int i=startRow;i<=endRow;i++) {
773 959145 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
774 319715 : uvw(2,i)=vb.uvw()(i)(2);
775 : }
776 463 : timemass_p +=tim.real();
777 463 : tim.mark();
778 463 : rotateUVW(uvw, dphase, vb);
779 463 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
780 : Bool delphase;
781 463 : Complex * phasorstor=phasor.getStorage(delphase);
782 463 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
783 463 : const Double * scalestor=uvScale.getStorage(del);
784 463 : const Double * offsetstor=uvOffset.getStorage(del);
785 463 : const Double* uvstor= uvw.getStorage(del);
786 463 : Int * locstor=loc.getStorage(del);
787 463 : Int * offstor=off.getStorage(del);
788 463 : const Double *dpstor=dphase.getStorage(del);
789 : //Vector<Double> f1=interpVisFreq_p.copy();
790 463 : Int nvchan=nvischan;
791 : Int irow;
792 463 : Double cinv=Double(1.0)/C::c;
793 463 : Int dow=0;
794 463 : Int nth=1;
795 : #ifdef _OPENMP
796 463 : if(numthreads_p >0){
797 0 : nth=min(numthreads_p, omp_get_max_threads());
798 : }
799 : else{
800 463 : nth= omp_get_max_threads();
801 : }
802 463 : nth=min(4,nth);
803 : #endif
804 :
805 :
806 463 : #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
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 : //////TESTOO
819 : //ofstream myfile;
820 : //myfile.open ("putLoc.txt", ios::out | ios::app | ios::ate );
821 : //myfile << vb.rowIds()(0) << " uv " << uvw.column(0) << " loc " << loc(0,0,0) << ", " << loc(1,0,0) << "\n" << endl;
822 : //myfile.close();
823 : ///////////////
824 :
825 463 : Int idopsf=0;
826 463 : if(dopsf) idopsf=1;
827 :
828 :
829 463 : Vector<Int> rowFlags(vb.nRow());
830 463 : rowFlags=0;
831 463 : rowFlags(vb.flagRow())=true;
832 463 : if(!usezero_p) {
833 320178 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
834 319715 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
835 : }
836 : }
837 :
838 : /////////////Some extra stuff for openmp
839 :
840 : Int x0, y0, nxsub, nysub, ixsub, iysub, icounter, ix, iy;
841 463 : Int csupp=convSupport_p;
842 :
843 463 : const Double * convfuncstor=(gridder->cFunction()).getStorage(del);
844 :
845 : // cerr <<"Poffset " << min(off) << " " << max(off) << " length " << gridder->cFunction().shape() << endl;
846 :
847 :
848 :
849 463 : ixsub=1;
850 463 : iysub=1;
851 463 : if (nth >3){
852 463 : ixsub=2;
853 463 : iysub=2;
854 : }
855 0 : else if(nth >1){
856 0 : ixsub=2;
857 0 : iysub=1;
858 : }
859 :
860 463 : x0=1;
861 463 : y0=1;
862 463 : nxsub=nx;
863 463 : nysub=ny;
864 463 : Int rbeg=startRow+1;
865 463 : Int rend=endRow+1;
866 463 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
867 2315 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
868 1852 : sumwgt[icounter].resize(sumWeight.shape());
869 1852 : sumwgt[icounter].set(0.0);
870 : }
871 463 : const Int* pmapstor=polMap.getStorage(del);
872 463 : const Int *cmapstor=chanMap.getStorage(del);
873 463 : Int nc=nchan;
874 463 : Int np=npol;
875 463 : Int nxp=nx;
876 463 : Int nyp=ny;
877 463 : const Int * flagstor=flags.getStorage(del);
878 463 : const Int * rowflagstor=rowFlags.getStorage(del);
879 : ////////////////////////
880 :
881 : Bool gridcopy;
882 463 : if(useDoubleGrid_p){
883 263 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
884 263 : #pragma omp parallel default(none) private(icounter,ix,iy,x0,y0,nxsub,nysub, 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) shared(sumwgt) num_threads(ixsub*iysub)
885 :
886 : {
887 : //cerr << "numthreads " << omp_get_num_threads() << endl;
888 : #pragma omp for
889 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
890 : //cerr << "thread id " << omp_get_thread_num() << endl;
891 : ix= (icounter+1)-((icounter)/ixsub)*ixsub;
892 : iy=(icounter)/ixsub+1;
893 : y0=(nyp/iysub)*(iy-1)+1;
894 : nysub=nyp/iysub;
895 : if( iy == iysub) {
896 : nysub=nyp-(nyp/iysub)*(iy-1);
897 : }
898 : x0=(nxp/ixsub)*(ix-1)+1;
899 : nxsub=nxp/ixsub;
900 : if( ix == ixsub){
901 : nxsub=nxp-(nxp/ixsub)*(ix-1);
902 : }
903 : sectggridd(datStorage,
904 : &nvispol,
905 : &nvischan,
906 : &idopsf,
907 : flagstor,
908 : rowflagstor,
909 : wgtStorage,
910 : &nvisrow,
911 : gridstor,
912 : &nxp,
913 : &nyp,
914 : &np,
915 : &nc,
916 : &csupp,
917 : &csamp,
918 : convfuncstor,
919 : cmapstor,
920 : pmapstor,
921 : (sumwgt[icounter]).getStorage(del),
922 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
923 : phasorstor);
924 : }//end for
925 : }// end pragma parallel
926 1315 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
927 1052 : sumWeight=sumWeight+sumwgt[icounter];
928 : }
929 : //phasor.putStorage(phasorstor, delphase);
930 263 : griddedData2.putStorage(gridstor, gridcopy);
931 : }
932 : else{
933 200 : Complex *gridstor=griddedData.getStorage(gridcopy);
934 200 : #pragma omp parallel default(none) private(icounter,ix,iy,x0,y0,nxsub,nysub, 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) shared(sumwgt) num_threads(ixsub*iysub)
935 : {
936 : //cerr << "numthreads " << omp_get_num_threads() << endl;
937 : #pragma omp for
938 :
939 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
940 : //cerr << "thread id " << omp_get_thread_num() << endl;
941 : ix= (icounter+1)-((icounter)/ixsub)*ixsub;
942 : iy=(icounter)/ixsub+1;
943 : y0=(nyp/iysub)*(iy-1)+1;
944 : nysub=nyp/iysub;
945 : if( iy == iysub) {
946 : nysub=nyp-(nyp/iysub)*(iy-1);
947 : }
948 : x0=(nxp/ixsub)*(ix-1)+1;
949 : nxsub=nxp/ixsub;
950 : if( ix == ixsub){
951 : nxsub=nxp-(nxp/ixsub)*(ix-1);
952 : }
953 : //cerr << "x0 " << x0 << " y0 " << y0 << " nxsub " << nxsub << " nysub " << nysub << endl;
954 : sectggrids(datStorage,
955 : &nvispol,
956 : &nvischan,
957 : &idopsf,
958 : flagstor,
959 : rowflagstor,
960 : wgtStorage,
961 : &nvisrow,
962 : gridstor,
963 : &nxp,
964 : &nyp,
965 : &np,
966 : &nc,
967 : &csupp,
968 : &csamp,
969 : convfuncstor,
970 : cmapstor,
971 : pmapstor,
972 : (sumwgt[icounter]).getStorage(del),
973 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
974 : phasorstor);
975 : }//end for
976 : }// end pragma parallel
977 1000 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
978 800 : sumWeight=sumWeight+sumwgt[icounter];
979 : }
980 :
981 200 : griddedData.putStorage(gridstor, gridcopy);
982 : }
983 : // cerr << "sunweight " << sumWeight << endl;
984 :
985 463 : timegrid_p+=tim.real();
986 :
987 463 : if(!dopsf)
988 200 : data.freeStorage(datStorage, isCopy);
989 463 : elWeight.freeStorage(wgtStorage,iswgtCopy);
990 :
991 : // peek->reset();
992 463 : }
993 :
994 0 : void GridFT::modifyConvFunc(const Vector<Double>& convFunc, Int convSupport, Int convSampling){
995 0 : convFunc_p.resize();
996 0 : convFunc_p=convFunc;
997 0 : convSupport_p=convSupport;
998 0 : convSampling_p=convSampling;
999 :
1000 0 : }
1001 :
1002 40 : void GridFT::get(VisBuffer& vb, Int row)
1003 : {
1004 :
1005 40 : gridOk(convSupport_p);
1006 : // If row is -1 then we pass through all rows
1007 : Int startRow, endRow, nRow;
1008 40 : if (row < 0) {
1009 40 : nRow=vb.nRow();
1010 40 : startRow=0;
1011 40 : endRow=nRow-1;
1012 : //unnecessary zeroing
1013 : //vb.modelVisCube()=Complex(0.0,0.0);
1014 : } else {
1015 0 : nRow=1;
1016 0 : startRow=row;
1017 0 : endRow=row;
1018 : //unnecessary zeroing
1019 : //vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
1020 : }
1021 :
1022 : // Get the uvws in a form that Fortran can use
1023 40 : Matrix<Double> uvw(3, vb.uvw().nelements());
1024 40 : uvw=0.0;
1025 40 : Vector<Double> dphase(vb.uvw().nelements());
1026 40 : dphase=0.0;
1027 : //NEGATING to correct for an image inversion problem
1028 40780 : for (Int i=startRow;i<=endRow;i++) {
1029 122220 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
1030 40740 : uvw(2,i)=vb.uvw()(i)(2);
1031 : }
1032 40 : rotateUVW(uvw, dphase, vb);
1033 40 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1034 :
1035 :
1036 :
1037 : //Check if ms has changed then cache new spw and chan selection
1038 40 : if(vb.newMS())
1039 0 : matchAllSpwChans(vb);
1040 :
1041 :
1042 : //Here we redo the match or use previous match
1043 :
1044 : //Channel matching for the actual spectral window of buffer
1045 40 : if(doConversion_p[vb.spectralWindow()]){
1046 0 : matchChannel(vb.spectralWindow(), vb);
1047 : }
1048 : else{
1049 40 : chanMap.resize();
1050 40 : chanMap=multiChanMap_p[vb.spectralWindow()];
1051 : }
1052 :
1053 : //cerr << "chanMap " << chanMap << endl;
1054 : //No point in reading data if its not matching in frequency
1055 40 : if(max(chanMap)==-1)
1056 0 : return;
1057 :
1058 40 : Cube<Complex> data;
1059 40 : Cube<Int> flags;
1060 40 : getInterpolateArrays(vb, data, flags);
1061 :
1062 : //cerr << "max min data " << max(griddedData) << " " << min(griddedData) << " flags " << max(flags) << " " << min(flags) << endl;
1063 :
1064 : // IPosition s(data.shape());
1065 40 : Int nvp=data.shape()(0);
1066 40 : Int nvc=data.shape()(1);
1067 40 : Int nvisrow=data.shape()(2);
1068 :
1069 :
1070 : //cerr << "get flags " << min(flags) << " " << max(flags) << endl;
1071 : Complex *datStorage;
1072 : Bool isCopy;
1073 40 : datStorage=data.getStorage(isCopy);
1074 :
1075 : ///
1076 40 : Cube<Int> loc(2, nvc, nvisrow);
1077 40 : Cube<Int> off(2, nvc, nRow);
1078 40 : Matrix<Complex> phasor(nvc, nRow);
1079 40 : Int csamp=convSampling_p;
1080 : Bool delphase;
1081 : Bool del;
1082 40 : Complex * phasorstor=phasor.getStorage(delphase);
1083 40 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1084 40 : const Double * scalestor=uvScale.getStorage(del);
1085 40 : const Double * offsetstor=uvOffset.getStorage(del);
1086 40 : const Double* uvstor= uvw.getStorage(del);
1087 40 : Int * locstor=loc.getStorage(del);
1088 40 : Int * offstor=off.getStorage(del);
1089 40 : const Double *dpstor=dphase.getStorage(del);
1090 : //Vector<Double> f1=interpVisFreq_p.copy();
1091 40 : Int nvchan=nvc;
1092 : Int irow;
1093 40 : Double cinv=Double(1.0)/C::c;
1094 40 : Int dow=0;
1095 40 : Int nth=1;
1096 : #ifdef _OPENMP
1097 40 : if(numthreads_p >0){
1098 0 : nth=min(numthreads_p, omp_get_max_threads());
1099 : }
1100 : else{
1101 40 : nth= omp_get_max_threads();
1102 : }
1103 40 : nth=min(4,nth);
1104 : #endif
1105 :
1106 :
1107 40 : #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)
1108 :
1109 : {
1110 : #pragma omp for
1111 : for (irow=startRow; irow<=endRow; ++irow){
1112 : //locateuvw(uvstor,dpstor, visfreqstor, nvchan, scalestor, offsetstor, csamp,
1113 : // locstor,
1114 : // offstor, phasorstor, irow);
1115 :
1116 : locuvw(uvstor, dpstor, visfreqstor, &nvchan, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
1117 : }
1118 :
1119 : }//end pragma parallel
1120 :
1121 40 : Int rbeg=startRow+1;
1122 40 : Int rend=endRow+1;
1123 :
1124 :
1125 40 : Vector<Int> rowFlags(vb.nRow());
1126 40 : rowFlags=0;
1127 40 : rowFlags(vb.flagRow())=true;
1128 : //cerr << "rowFlags " << rowFlags << endl;
1129 40 : if(!usezero_p) {
1130 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1131 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
1132 : }
1133 : }
1134 :
1135 :
1136 : //cerr <<"offset " << min(off) << " " <<max(off) << " length " << gridder->cFunction().shape() << endl;
1137 :
1138 : {
1139 : Bool delgrid;
1140 40 : const Complex* gridstor=griddedData.getStorage(delgrid);
1141 40 : const Double * convfuncstor=convFunc_p.getStorage(del);;
1142 :
1143 40 : const Int* pmapstor=polMap.getStorage(del);
1144 40 : const Int *cmapstor=chanMap.getStorage(del);
1145 40 : Int nc=nchan;
1146 40 : Int np=npol;
1147 40 : Int nxp=nx;
1148 40 : Int nyp=ny;
1149 40 : Int csupp=convSupport_p;
1150 40 : const Int * flagstor=flags.getStorage(del);
1151 40 : const Int * rowflagstor=rowFlags.getStorage(del);
1152 :
1153 :
1154 40 : Int npart=1;
1155 40 : if (nth >3){
1156 40 : npart=4;
1157 : }
1158 0 : else if(nth >1){
1159 0 : npart=2;
1160 : }
1161 :
1162 40 : Int ix=0;
1163 40 : #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)
1164 : {
1165 : #pragma omp for
1166 : for (ix=0; ix< npart; ++ix){
1167 : rbeg=ix*(nvisrow/npart)+1;
1168 : rend=(ix != (npart-1)) ? (rbeg+(nvisrow/npart)-1) : (rbeg+(nvisrow/npart)+nvisrow%npart-1) ;
1169 : //cerr << "rbeg " << rbeg << " rend " << rend << " " << nvisrow << endl;
1170 :
1171 : sectdgrid(datStorage,
1172 : &nvp,
1173 : &nvc,
1174 : flagstor,
1175 : rowflagstor,
1176 : &nvisrow,
1177 : gridstor,
1178 : &nxp,
1179 : &nyp,
1180 : &np,
1181 : &nc,
1182 : &csupp,
1183 : &csamp,
1184 : convfuncstor,
1185 : cmapstor,
1186 : pmapstor,
1187 : &rbeg, &rend, locstor, offstor, phasorstor);
1188 : }//end pragma parallel
1189 : }
1190 40 : data.putStorage(datStorage, isCopy);
1191 40 : griddedData.freeStorage(gridstor, delgrid);
1192 : //cerr << "Get min max " << min(data) << " " << max(data) << endl;
1193 : }
1194 40 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1195 :
1196 40 : }
1197 :
1198 :
1199 :
1200 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1201 : // return the resulting image
1202 15 : ImageInterface<Complex>& GridFT::getImage(Matrix<Float>& weights, Bool normalize)
1203 : {
1204 : //AlwaysAssert(lattice, AipsError);
1205 15 : AlwaysAssert(gridder, AipsError);
1206 15 : AlwaysAssert(image, AipsError);
1207 15 : logIO() << LogOrigin("GridFT", "getImage") << LogIO::NORMAL;
1208 :
1209 15 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
1210 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
1211 : }
1212 15 : weights.resize(sumWeight.shape());
1213 :
1214 15 : convertArray(weights, sumWeight);
1215 : // If the weights are all zero then we cannot normalize
1216 : // otherwise we don't care.
1217 15 : if(normalize&&max(weights)==0.0) {
1218 : logIO() << LogIO::SEVERE << "No useful data in GridFT: weights all zero"
1219 0 : << LogIO::POST;
1220 : }
1221 : else {
1222 :
1223 : //const IPosition latticeShape = lattice->shape();
1224 :
1225 : //logIO() << LogIO::DEBUGGING
1226 : // << "Starting FFT and scaling of image" << LogIO::POST;
1227 :
1228 :
1229 :
1230 : // if(useDoubleGrid_p){
1231 : // convertArray(griddedData, griddedData2);
1232 : // //Don't need the double-prec grid anymore...
1233 : // griddedData2.resize();
1234 : // }
1235 :
1236 : // x and y transforms
1237 : // LatticeFFT::cfft2d(*lattice,false);
1238 : //
1239 : // Retain the double precision grid for FFT as well. Convert it
1240 : // to single precision just after (since images are still single
1241 : // precision).
1242 : //
1243 15 : if(useDoubleGrid_p)
1244 : {
1245 14 : ArrayLattice<DComplex> darrayLattice(griddedData2);
1246 :
1247 14 : LatticeFFT::cfft2d(darrayLattice,false);
1248 14 : griddedData.resize(griddedData2.shape());
1249 14 : convertArray(griddedData, griddedData2);
1250 :
1251 14 : SynthesisUtilMethods::getResource("mem peak in getImage");
1252 :
1253 : //Don't need the double-prec grid anymore...
1254 14 : griddedData2.resize();
1255 14 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1256 14 : lattice=arrayLattice;
1257 14 : }
1258 : else{
1259 1 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1260 : //storeArrayAsImage(String("cgrid_ft.im"), image->coordinates(), griddedData);
1261 1 : lattice=arrayLattice;
1262 1 : LatticeFFT::cfft2d(*lattice,false);
1263 : }
1264 :
1265 :
1266 :
1267 : {
1268 15 : Int inx = lattice->shape()(0);
1269 15 : Int iny = lattice->shape()(1);
1270 15 : Vector<Complex> correction(inx);
1271 15 : correction=Complex(1.0, 0.0);
1272 : // Do the Grid-correction
1273 15 : IPosition cursorShape(4, inx, 1, 1, 1);
1274 15 : IPosition axisPath(4, 0, 1, 2, 3);
1275 15 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1276 15 : LatticeIterator<Complex> lix(*lattice, lsx);
1277 2835 : for(lix.reset();!lix.atEnd();lix++) {
1278 2820 : Int pol=lix.position()(2);
1279 2820 : Int chan=lix.position()(3);
1280 2820 : if(weights(pol, chan)!=0.0) {
1281 2820 : gridder->correctX1D(correction, lix.position()(1));
1282 2820 : lix.rwVectorCursor()/=correction;
1283 2820 : if(normalize) {
1284 100 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
1285 100 : lix.rwCursor()*=rnorm;
1286 : }
1287 : else {
1288 2720 : Complex rnorm(Float(inx)*Float(iny));
1289 2720 : lix.rwCursor()*=rnorm;
1290 : }
1291 : }
1292 : else {
1293 0 : lix.woCursor()=0.0;
1294 : }
1295 : }
1296 15 : }
1297 :
1298 15 : if(!isTiled) {
1299 : // Check the section from the image BEFORE converting to a lattice
1300 30 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1301 15 : IPosition stride(4, 1);
1302 30 : IPosition trc(blc+image->shape()-stride);
1303 : // Do the copy
1304 15 : IPosition start(4, 0);
1305 15 : image->put(griddedData(blc, trc));
1306 15 : }
1307 : }
1308 :
1309 15 : image->flush();
1310 15 : image->clearCache();
1311 :
1312 15 : if(!arrayLattice.null()) arrayLattice=0;
1313 15 : if(!lattice.null()) lattice=0;
1314 15 : griddedData.resize();
1315 :
1316 15 : return *image;
1317 : }
1318 :
1319 : // Get weight image
1320 14 : void GridFT::getWeightImage(ImageInterface<Float>& weightImage, Matrix<Float>& weights)
1321 : {
1322 :
1323 14 : logIO() << LogOrigin("GridFT", "getWeightImage") << LogIO::NORMAL;
1324 :
1325 14 : weights.resize(sumWeight.shape());
1326 14 : convertArray(weights,sumWeight);
1327 :
1328 14 : const IPosition latticeShape = weightImage.shape();
1329 :
1330 14 : Int nx=latticeShape(0);
1331 14 : Int ny=latticeShape(1);
1332 :
1333 14 : IPosition loc(2, 0);
1334 14 : IPosition cursorShape(4, nx, ny, 1, 1);
1335 14 : IPosition axisPath(4, 0, 1, 2, 3);
1336 14 : LatticeStepper lsx(latticeShape, cursorShape, axisPath);
1337 14 : LatticeIterator<Float> lix(weightImage, lsx);
1338 28 : for(lix.reset();!lix.atEnd();lix++) {
1339 14 : Int pol=lix.position()(2);
1340 14 : Int chan=lix.position()(3);
1341 14 : lix.rwCursor()=weights(pol,chan);
1342 : }
1343 14 : }
1344 :
1345 0 : Bool GridFT::toRecord(String& error,
1346 : RecordInterface& outRec, Bool withImage, const String diskimage)
1347 : {
1348 :
1349 :
1350 :
1351 : // Save the current GridFT object to an output state record
1352 0 : Bool retval = true;
1353 0 : Float elpadd=padding_p;
1354 0 : if(toVis_p && withImage)
1355 0 : elpadd=1.0;
1356 : //save the base class variables
1357 0 : if(!FTMachine::toRecord(error, outRec, withImage, diskimage))
1358 0 : return false;
1359 :
1360 : //a call to init will redo imagecache and gridder
1361 : // so no need to save these
1362 :
1363 0 : outRec.define("cachesize", Int64(cachesize));
1364 0 : outRec.define("tilesize", tilesize);
1365 0 : outRec.define("convtype", convType);
1366 0 : outRec.define("uvscale", uvScale);
1367 0 : outRec.define("uvoffset", uvOffset);
1368 0 : outRec.define("istiled", isTiled);
1369 0 : outRec.define("maxabsdata", maxAbsData);
1370 0 : Vector<Int> center_loc(4), offset_loc(4);
1371 0 : for (Int k=0; k<4 ; k++){
1372 0 : center_loc(k)=centerLoc(k);
1373 0 : offset_loc(k)=offsetLoc(k);
1374 : }
1375 0 : outRec.define("centerloc", center_loc);
1376 0 : outRec.define("offsetloc", offset_loc);
1377 0 : outRec.define("padding", elpadd);
1378 0 : outRec.define("usezero", usezero_p);
1379 0 : outRec.define("nopadding", noPadding_p);
1380 0 : return retval;
1381 0 : }
1382 :
1383 0 : Bool GridFT::fromRecord(String& error,
1384 : const RecordInterface& inRec)
1385 : {
1386 0 : Bool retval = true;
1387 0 : if(!FTMachine::fromRecord(error, inRec))
1388 0 : return false;
1389 0 : gridder=0; imageCache=0; lattice=0; arrayLattice=0;
1390 : //For some reason int64 seems to be getting lost...
1391 : //this is not an important parameter...so filing a jira
1392 0 : if(inRec.isDefined("cachesize"))
1393 0 : cachesize=inRec.asInt64("cachesize");
1394 : else
1395 0 : cachesize=1000000;
1396 0 : inRec.get("tilesize", tilesize);
1397 0 : inRec.get("convtype", convType);
1398 0 : inRec.get("uvscale", uvScale);
1399 0 : inRec.get("uvoffset", uvOffset);
1400 0 : inRec.get("istiled", isTiled);
1401 0 : inRec.get("maxabsdata", maxAbsData);
1402 0 : Vector<Int> center_loc(4), offset_loc(4);
1403 0 : inRec.get("centerloc", center_loc);
1404 0 : inRec.get("offsetloc", offset_loc);
1405 0 : uInt ndim4 = 4;
1406 0 : centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2),
1407 0 : center_loc(3));
1408 0 : offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2),
1409 0 : offset_loc(3));
1410 0 : inRec.get("padding", padding_p);
1411 0 : inRec.get("usezero", usezero_p);
1412 0 : inRec.get("nopadding", noPadding_p);
1413 :
1414 0 : machineName_p="GridFT";
1415 : ///setup some of the parameters
1416 0 : init();
1417 : /*if(!cmplxImage_p.null()){
1418 : //FTMachine::fromRecord would have recovered the image
1419 : // Might be changing the shape of sumWeight
1420 :
1421 : if(isTiled) {
1422 : lattice=CountedPtr<Lattice<Complex> >(image, false);
1423 : }
1424 : else {
1425 : // Make the grid the correct shape and turn it into an array lattice
1426 : // Check the section from the image BEFORE converting to a lattice
1427 : if(!toVis_p){
1428 : IPosition gridShape(4, nx, ny, npol, nchan);
1429 : griddedData.resize(gridShape);
1430 : griddedData=Complex(0.0);
1431 : }
1432 : else{
1433 : prepGridForDegrid();
1434 : }
1435 : }
1436 : };*/
1437 0 : return retval;
1438 0 : }
1439 :
1440 1101 : void GridFT::ok() {
1441 1101 : AlwaysAssert(image, AipsError);
1442 1101 : }
1443 :
1444 : // Make a plain straightforward honest-to-God image. This returns
1445 : // a complex image, without conversion to Stokes. The representation
1446 : // is that required for the visibilities.
1447 : //----------------------------------------------------------------------
1448 0 : void GridFT::makeImage(FTMachine::Type type,
1449 : VisSet& vs,
1450 : ImageInterface<Complex>& theImage,
1451 : Matrix<Float>& weight) {
1452 :
1453 :
1454 0 : logIO() << LogOrigin("GridFT", "makeImage") << LogIO::NORMAL;
1455 :
1456 0 : if(type==FTMachine::COVERAGE) {
1457 0 : logIO() << "Type COVERAGE not defined for Fourier transforms" << LogIO::EXCEPTION;
1458 : }
1459 :
1460 :
1461 : // Initialize the gradients
1462 0 : ROVisIter& vi(vs.iter());
1463 :
1464 : // Loop over all visibilities and pixels
1465 0 : VisBuffer vb(vi);
1466 :
1467 : // Initialize put (i.e. transform to Sky) for this model
1468 0 : vi.origin();
1469 :
1470 0 : if(vb.polFrame()==MSIter::Linear) {
1471 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
1472 : }
1473 : else {
1474 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
1475 : }
1476 :
1477 0 : initializeToSky(theImage,weight,vb);
1478 :
1479 : // Loop over the visibilities, putting VisBuffers
1480 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1481 0 : for (vi.origin(); vi.more(); vi++) {
1482 :
1483 0 : switch(type) {
1484 0 : case FTMachine::RESIDUAL:
1485 0 : vb.visCube()=vb.correctedVisCube();
1486 0 : vb.visCube()-=vb.modelVisCube();
1487 0 : put(vb, -1, false);
1488 0 : break;
1489 0 : case FTMachine::MODEL:
1490 0 : vb.visCube()=vb.modelVisCube();
1491 0 : put(vb, -1, false);
1492 0 : break;
1493 0 : case FTMachine::CORRECTED:
1494 0 : vb.visCube()=vb.correctedVisCube();
1495 0 : put(vb, -1, false);
1496 0 : break;
1497 0 : case FTMachine::PSF:
1498 0 : vb.visCube()=Complex(1.0,0.0);
1499 0 : put(vb, -1, true);
1500 0 : break;
1501 0 : case FTMachine::OBSERVED:
1502 : default:
1503 0 : put(vb, -1, false);
1504 0 : break;
1505 : }
1506 : }
1507 : }
1508 0 : finalizeToSky();
1509 : // Normalize by dividing out weights, etc.
1510 0 : getImage(weight, true);
1511 0 : }
1512 :
1513 3740 : String GridFT::name() const {
1514 :
1515 3740 : return machineName_p;
1516 :
1517 :
1518 : }
1519 :
1520 : } //# NAMESPACE CASA - END
1521 :
|
