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: aips2-request@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 8 : GridFT::GridFT() : FTMachine(), padding_p(1.0), imageCache(0), cachesize(1000000), tilesize(1000), gridder(0), isTiled(false), convType("SF"),
83 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
84 : usezero_p(false), noPadding_p(false), usePut2_p(false),
85 8 : machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0) {
86 :
87 8 : }
88 0 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType, Float padding,
89 0 : Bool usezero, Bool useDoublePrec)
90 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize), tilesize(itilesize),
91 : gridder(0), isTiled(false), convType(iconvType),
92 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
93 : 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 143 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
101 143 : MPosition mLocation, Float padding, Bool usezero, Bool useDoublePrec)
102 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
103 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
104 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
105 143 : 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 143 : mLocation_p=mLocation;
108 143 : tangentSpecified_p=false;
109 143 : useDoubleGrid_p=useDoublePrec;
110 : // peek=NULL;
111 143 : }
112 :
113 0 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
114 0 : MDirection mTangent, Float padding, Bool usezero, Bool useDoublePrec)
115 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
116 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
117 : 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 4 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
127 : MPosition mLocation, MDirection mTangent, Float padding,
128 4 : Bool usezero, Bool useDoublePrec)
129 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
130 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
131 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
132 4 : 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 4 : mLocation_p=mLocation;
135 4 : mTangent_p=mTangent;
136 4 : tangentSpecified_p=true;
137 4 : useDoubleGrid_p=useDoublePrec;
138 : // peek=NULL;
139 4 : }
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 262 : GridFT& GridFT::operator=(const GridFT& other)
157 : {
158 262 : if(this!=&other) {
159 : //Do the base parameters
160 262 : FTMachine::operator=(other);
161 :
162 : //private params
163 262 : imageCache=other.imageCache;
164 262 : cachesize=other.cachesize;
165 262 : tilesize=other.tilesize;
166 262 : convType=other.convType;
167 262 : convType.upcase();
168 262 : uvScale.resize();
169 262 : uvOffset.resize();
170 262 : uvScale.assign(other.uvScale);
171 262 : uvOffset.assign(other.uvOffset);
172 262 : if(other.gridder==0)
173 262 : gridder=0;
174 : else{
175 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
176 0 : uvScale, uvOffset,
177 0 : convType);
178 : }
179 262 : convFunc_p.resize();
180 262 : convSupport_p=other.convSupport_p;
181 262 : convSampling_p=other.convSampling_p;
182 262 : convFunc_p=other.convFunc_p;
183 262 : isTiled=other.isTiled;
184 : //lattice=other.lattice;
185 262 : lattice=0;
186 262 : tilesize=other.tilesize;
187 262 : arrayLattice=0;
188 262 : maxAbsData=other.maxAbsData;
189 262 : centerLoc=other.centerLoc;
190 262 : offsetLoc=other.offsetLoc;
191 262 : padding_p=other.padding_p;
192 262 : usezero_p=other.usezero_p;
193 262 : noPadding_p=other.noPadding_p;
194 262 : machineName_p="GridFT";
195 262 : timemass_p=0.0;
196 262 : timegrid_p=0.0;
197 : // peek = other.peek;
198 : };
199 262 : return *this;
200 : };
201 :
202 : //----------------------------------------------------------------------
203 254 : GridFT::GridFT(const GridFT& other) : FTMachine(), machineName_p("GridFT")
204 : {
205 254 : operator=(other);
206 254 : }
207 :
208 : //-----------------------------------------------------------------------
209 0 : FTMachine* GridFT::cloneFTM(){
210 0 : return new GridFT(*this);
211 : }
212 :
213 : //----------------------------------------------------------------------
214 69 : void GridFT::init() {
215 :
216 69 : logIO() << LogOrigin("GridFT", "init") << LogIO::NORMAL;
217 :
218 69 : 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 69 : isTiled=false;
239 69 : if(!noPadding_p){
240 69 : CompositeNumber cn(uInt(image->shape()(0)*2));
241 69 : nx = cn.nextLargerEven(Int(padding_p*Float(image->shape()(0))-0.5));
242 69 : ny = cn.nextLargerEven(Int(padding_p*Float(image->shape()(1))-0.5));
243 : }
244 : else{
245 0 : nx = image->shape()(0);
246 0 : ny = image->shape()(1);
247 : }
248 69 : npol = image->shape()(2);
249 69 : nchan = image->shape()(3);
250 : // }
251 :
252 69 : sumWeight.resize(npol, nchan);
253 :
254 69 : uvScale.resize(2);
255 69 : uvScale(0)=(Float(nx)*image->coordinates().increment()(0));
256 69 : uvScale(1)=(Float(ny)*image->coordinates().increment()(1));
257 69 : uvOffset.resize(2);
258 69 : uvOffset(0)=nx/2;
259 69 : uvOffset(1)=ny/2;
260 :
261 : // Now set up the gridder. The possibilities are BOX and SF
262 69 : if(gridder) delete gridder; gridder=0;
263 69 : convType.upcase();
264 138 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
265 69 : uvScale, uvOffset,
266 69 : convType);
267 :
268 :
269 69 : convFunc_p.resize();
270 69 : convSupport_p=gridder->cSupport()(0);
271 69 : convSampling_p=gridder->cSampling();
272 69 : 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 69 : if(imageCache) delete imageCache; imageCache=0;
278 :
279 69 : 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 : }
299 69 : }
300 :
301 : // This is nasty, we should use CountedPointers here.
302 806 : GridFT::~GridFT() {
303 409 : if(imageCache) delete imageCache; imageCache=0;
304 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
305 409 : if(gridder) delete gridder; gridder=0;
306 806 : }
307 :
308 : // Initialize for a transform from the Sky domain. This means that
309 : // we grid-correct, and FFT the image
310 :
311 55 : void GridFT::initializeToVis(ImageInterface<Complex>& iimage,
312 : const VisBuffer& vb)
313 : {
314 55 : image=&iimage;
315 55 : toVis_p=true;
316 :
317 55 : ok();
318 :
319 55 : init();
320 : // peek->reset();
321 : // Initialize the maps for polarization and channel. These maps
322 : // translate visibility indices into image indices
323 55 : 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 55 : prepGridForDegrid();
342 :
343 : //AlwaysAssert(lattice, AipsError);
344 :
345 :
346 :
347 55 : }
348 :
349 55 : void GridFT::prepGridForDegrid(){
350 110 : IPosition gridShape(4, nx, ny, npol, nchan);
351 55 : 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 55 : griddedData.set(Complex(0.0));
357 :
358 110 : IPosition stride(4, 1);
359 165 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
360 165 : IPosition trc(blc+image->shape()-stride);
361 :
362 110 : IPosition start(4, 0);
363 55 : griddedData(blc, trc) = image->getSlice(start, image->shape());
364 :
365 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
366 55 : arrayLattice = new ArrayLattice<Complex>(griddedData);
367 55 : lattice=arrayLattice;
368 : //logIO() << LogIO::DEBUGGING
369 : // << "Starting grid correction and FFT of image" << LogIO::POST;
370 :
371 : // Do the Grid-correction.
372 : {
373 110 : Vector<Complex> correction(nx);
374 55 : correction=Complex(1.0, 0.0);
375 : // Do the Grid-correction
376 110 : IPosition cursorShape(4, nx, 1, 1, 1);
377 110 : IPosition axisPath(4, 0, 1, 2, 3);
378 110 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
379 110 : LatticeIterator<Complex> lix(*lattice, lsx);
380 112551 : for(lix.reset();!lix.atEnd();lix++) {
381 112496 : gridder->correctX1D(correction, lix.position()(1));
382 112496 : lix.rwVectorCursor()/=correction;
383 : }
384 : }
385 55 : image->clearCache();
386 : // Now do the FFT2D in place
387 55 : LatticeFFT::cfft2d(*lattice);
388 :
389 : //logIO() << LogIO::DEBUGGING
390 : // << "Finished grid correction and FFT of image" << LogIO::POST;
391 :
392 55 : }
393 :
394 :
395 12 : void GridFT::finalizeToVis()
396 : {
397 :
398 : //cerr <<"Time to degrid " << timedegrid_p << endl;
399 12 : timedegrid_p=0.0;
400 :
401 12 : if(!arrayLattice.null()) arrayLattice=0;
402 12 : if(!lattice.null()) lattice=0;
403 12 : griddedData.resize();
404 :
405 12 : 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 : }
416 12 : }
417 :
418 :
419 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
420 : // grid.
421 14 : void GridFT::initializeToSky(ImageInterface<Complex>& iimage,
422 : Matrix<Float>& weight, const VisBuffer& vb)
423 : {
424 : // image always points to the image
425 14 : image=&iimage;
426 14 : toVis_p=false;
427 :
428 14 : init();
429 :
430 : // Initialize the maps for polarization and channel. These maps
431 : // translate visibility indices into image indices
432 14 : initMaps(vb);
433 :
434 :
435 :
436 14 : sumWeight=0.0;
437 14 : weight.resize(sumWeight.shape());
438 14 : 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 14 : 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 28 : IPosition gridShape(4, nx, ny, npol, nchan);
450 14 : if( !useDoubleGrid_p )
451 : {
452 0 : griddedData.resize(gridShape);
453 0 : griddedData=Complex(0.0);
454 : }
455 14 : if(useDoubleGrid_p){
456 : //griddedData.resize();
457 14 : griddedData2.resize(gridShape);
458 14 : griddedData2=DComplex(0.0);
459 : }
460 14 : image->clearCache();
461 : //iimage.get(griddedData, false);
462 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
463 :
464 : }
465 : //AlwaysAssert(lattice, AipsError);
466 14 : }
467 :
468 :
469 :
470 14 : 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 14 : timemass_p=0.0;
478 14 : timegrid_p=0.0;
479 14 : 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 : }
489 : // peek->terminate();
490 : // peek->reset();
491 14 : }
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 263 : void GridFT::put(const VisBuffer& vb, Int row, Bool dopsf,
682 : FTMachine::Type type)
683 : {
684 :
685 263 : gridOk(convSupport_p);
686 : // peek->setVBPtr(&vb);
687 :
688 : //Check if ms has changed then cache new spw and chan selection
689 263 : 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 263 : if(doConversion_p[vb.spectralWindow()]){
696 0 : matchChannel(vb.spectralWindow(), vb);
697 : }
698 : else{
699 263 : chanMap.resize();
700 263 : chanMap=multiChanMap_p[vb.spectralWindow()];
701 : }
702 :
703 : //No point in reading data if its not matching in frequency
704 263 : if(max(chanMap)==-1)
705 : {
706 : // peek->reset();
707 0 : return;
708 : }
709 :
710 263 : Timer tim;
711 263 : tim.mark();
712 :
713 : const Matrix<Float> *imagingweight;
714 263 : imagingweight=&(vb.imagingWeight());
715 :
716 263 : if(dopsf) {type=FTMachine::PSF;}
717 :
718 526 : Cube<Complex> data;
719 : //Fortran gridder need the flag as ints
720 526 : Cube<Int> flags;
721 526 : Matrix<Float> elWeight;
722 263 : interpolateFrequencyTogrid(vb, *imagingweight,data, flags, elWeight, type);
723 :
724 :
725 : Bool iswgtCopy;
726 : const Float *wgtStorage;
727 263 : wgtStorage=elWeight.getStorage(iswgtCopy);
728 :
729 : Bool isCopy;
730 : const Complex *datStorage;
731 263 : if(!dopsf)
732 0 : 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 263 : if (row==-1) {
738 263 : nRow=vb.nRow();
739 263 : startRow=0;
740 263 : 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 263 : const IPosition &fs=flags.shape();
759 526 : std::vector<Int> s(fs.begin(),fs.end());
760 263 : Int nvispol=s[0];
761 263 : Int nvischan=s[1];
762 263 : Int nvisrow=s[2];
763 526 : Matrix<Double> uvw(3, vb.uvw().nelements());
764 263 : uvw=0.0;
765 526 : Vector<Double> dphase(vb.uvw().nelements());
766 526 : Cube<Int> loc(2, nvischan, nRow);
767 526 : Cube<Int> off(2, nvischan, nRow);
768 526 : Matrix<Complex> phasor(nvischan, nRow);
769 263 : Int csamp=convSampling_p;
770 263 : dphase=0.0;
771 : //NEGATING to correct for an image inversion problem
772 249778 : for (Int i=startRow;i<=endRow;i++) {
773 748545 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
774 249515 : uvw(2,i)=vb.uvw()(i)(2);
775 : }
776 263 : timemass_p +=tim.real();
777 263 : tim.mark();
778 263 : rotateUVW(uvw, dphase, vb);
779 263 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
780 : Bool delphase;
781 263 : Complex * phasorstor=phasor.getStorage(delphase);
782 263 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
783 263 : const Double * scalestor=uvScale.getStorage(del);
784 263 : const Double * offsetstor=uvOffset.getStorage(del);
785 263 : const Double* uvstor= uvw.getStorage(del);
786 263 : Int * locstor=loc.getStorage(del);
787 263 : Int * offstor=off.getStorage(del);
788 263 : const Double *dpstor=dphase.getStorage(del);
789 : //Vector<Double> f1=interpVisFreq_p.copy();
790 263 : Int nvchan=nvischan;
791 : Int irow;
792 263 : Double cinv=Double(1.0)/C::c;
793 263 : Int dow=0;
794 263 : Int nth=1;
795 : #ifdef _OPENMP
796 263 : if(numthreads_p >0){
797 0 : nth=min(numthreads_p, omp_get_max_threads());
798 : }
799 : else{
800 263 : nth= omp_get_max_threads();
801 : }
802 263 : nth=min(4,nth);
803 : #endif
804 :
805 :
806 263 : #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 263 : Int idopsf=0;
826 263 : if(dopsf) idopsf=1;
827 :
828 :
829 526 : Vector<Int> rowFlags(vb.nRow());
830 263 : rowFlags=0;
831 263 : rowFlags(vb.flagRow())=true;
832 263 : if(!usezero_p) {
833 249778 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
834 249515 : 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 263 : Int csupp=convSupport_p;
842 :
843 263 : const Double * convfuncstor=(gridder->cFunction()).getStorage(del);
844 :
845 : // cerr <<"Poffset " << min(off) << " " << max(off) << " length " << gridder->cFunction().shape() << endl;
846 :
847 :
848 :
849 263 : ixsub=1;
850 263 : iysub=1;
851 263 : if (nth >3){
852 0 : ixsub=2;
853 0 : iysub=2;
854 : }
855 263 : else if(nth >1){
856 0 : ixsub=2;
857 0 : iysub=1;
858 : }
859 :
860 263 : x0=1;
861 263 : y0=1;
862 263 : nxsub=nx;
863 263 : nysub=ny;
864 263 : Int rbeg=startRow+1;
865 263 : Int rend=endRow+1;
866 526 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
867 526 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
868 263 : sumwgt[icounter].resize(sumWeight.shape());
869 263 : sumwgt[icounter].set(0.0);
870 : }
871 263 : const Int* pmapstor=polMap.getStorage(del);
872 263 : const Int *cmapstor=chanMap.getStorage(del);
873 263 : Int nc=nchan;
874 263 : Int np=npol;
875 263 : Int nxp=nx;
876 263 : Int nyp=ny;
877 263 : const Int * flagstor=flags.getStorage(del);
878 263 : const Int * rowflagstor=rowFlags.getStorage(del);
879 : ////////////////////////
880 :
881 : Bool gridcopy;
882 263 : 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 526 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
927 263 : sumWeight=sumWeight+sumwgt[icounter];
928 : }
929 : //phasor.putStorage(phasorstor, delphase);
930 263 : griddedData2.putStorage(gridstor, gridcopy);
931 : }
932 : else{
933 0 : Complex *gridstor=griddedData.getStorage(gridcopy);
934 0 : #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 0 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
978 0 : sumWeight=sumWeight+sumwgt[icounter];
979 : }
980 :
981 0 : griddedData.putStorage(gridstor, gridcopy);
982 : }
983 : // cerr << "sunweight " << sumWeight << endl;
984 :
985 263 : timegrid_p+=tim.real();
986 :
987 263 : if(!dopsf)
988 0 : data.freeStorage(datStorage, isCopy);
989 263 : elWeight.freeStorage(wgtStorage,iswgtCopy);
990 :
991 : // peek->reset();
992 : }
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 266 : void GridFT::get(VisBuffer& vb, Int row)
1003 : {
1004 :
1005 266 : gridOk(convSupport_p);
1006 : // If row is -1 then we pass through all rows
1007 : Int startRow, endRow, nRow;
1008 266 : if (row < 0) {
1009 266 : nRow=vb.nRow();
1010 266 : startRow=0;
1011 266 : 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 266 : Matrix<Double> uvw(3, vb.uvw().nelements());
1024 266 : uvw=0.0;
1025 266 : Vector<Double> dphase(vb.uvw().nelements());
1026 266 : dphase=0.0;
1027 : //NEGATING to correct for an image inversion problem
1028 120224 : for (Int i=startRow;i<=endRow;i++) {
1029 359874 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
1030 119958 : uvw(2,i)=vb.uvw()(i)(2);
1031 : }
1032 266 : rotateUVW(uvw, dphase, vb);
1033 266 : 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 266 : 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 266 : if(doConversion_p[vb.spectralWindow()]){
1046 0 : matchChannel(vb.spectralWindow(), vb);
1047 : }
1048 : else{
1049 266 : chanMap.resize();
1050 266 : 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 266 : if(max(chanMap)==-1)
1056 0 : return;
1057 :
1058 532 : Cube<Complex> data;
1059 532 : Cube<Int> flags;
1060 266 : 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 266 : Int nvp=data.shape()(0);
1066 266 : Int nvc=data.shape()(1);
1067 266 : Int nvisrow=data.shape()(2);
1068 :
1069 :
1070 : //cerr << "get flags " << min(flags) << " " << max(flags) << endl;
1071 : Complex *datStorage;
1072 : Bool isCopy;
1073 266 : datStorage=data.getStorage(isCopy);
1074 :
1075 : ///
1076 532 : Cube<Int> loc(2, nvc, nvisrow);
1077 532 : Cube<Int> off(2, nvc, nRow);
1078 532 : Matrix<Complex> phasor(nvc, nRow);
1079 266 : Int csamp=convSampling_p;
1080 : Bool delphase;
1081 : Bool del;
1082 266 : Complex * phasorstor=phasor.getStorage(delphase);
1083 266 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1084 266 : const Double * scalestor=uvScale.getStorage(del);
1085 266 : const Double * offsetstor=uvOffset.getStorage(del);
1086 266 : const Double* uvstor= uvw.getStorage(del);
1087 266 : Int * locstor=loc.getStorage(del);
1088 266 : Int * offstor=off.getStorage(del);
1089 266 : const Double *dpstor=dphase.getStorage(del);
1090 : //Vector<Double> f1=interpVisFreq_p.copy();
1091 266 : Int nvchan=nvc;
1092 : Int irow;
1093 266 : Double cinv=Double(1.0)/C::c;
1094 266 : Int dow=0;
1095 266 : Int nth=1;
1096 : #ifdef _OPENMP
1097 266 : if(numthreads_p >0){
1098 0 : nth=min(numthreads_p, omp_get_max_threads());
1099 : }
1100 : else{
1101 266 : nth= omp_get_max_threads();
1102 : }
1103 266 : nth=min(4,nth);
1104 : #endif
1105 :
1106 :
1107 266 : #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 266 : Int rbeg=startRow+1;
1122 266 : Int rend=endRow+1;
1123 :
1124 :
1125 532 : Vector<Int> rowFlags(vb.nRow());
1126 266 : rowFlags=0;
1127 266 : rowFlags(vb.flagRow())=true;
1128 : //cerr << "rowFlags " << rowFlags << endl;
1129 266 : if(!usezero_p) {
1130 68884 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1131 68688 : 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 266 : const Complex* gridstor=griddedData.getStorage(delgrid);
1141 266 : const Double * convfuncstor=convFunc_p.getStorage(del);;
1142 :
1143 266 : const Int* pmapstor=polMap.getStorage(del);
1144 266 : const Int *cmapstor=chanMap.getStorage(del);
1145 266 : Int nc=nchan;
1146 266 : Int np=npol;
1147 266 : Int nxp=nx;
1148 266 : Int nyp=ny;
1149 266 : Int csupp=convSupport_p;
1150 266 : const Int * flagstor=flags.getStorage(del);
1151 266 : const Int * rowflagstor=rowFlags.getStorage(del);
1152 :
1153 :
1154 266 : Int npart=1;
1155 266 : if (nth >3){
1156 0 : npart=4;
1157 : }
1158 266 : else if(nth >1){
1159 0 : npart=2;
1160 : }
1161 :
1162 266 : Int ix=0;
1163 266 : #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 266 : data.putStorage(datStorage, isCopy);
1191 266 : griddedData.freeStorage(gridstor, delgrid);
1192 : //cerr << "Get min max " << min(data) << " " << max(data) << endl;
1193 : }
1194 266 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1195 :
1196 : }
1197 :
1198 :
1199 :
1200 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1201 : // return the resulting image
1202 14 : ImageInterface<Complex>& GridFT::getImage(Matrix<Float>& weights, Bool normalize)
1203 : {
1204 : //AlwaysAssert(lattice, AipsError);
1205 14 : AlwaysAssert(gridder, AipsError);
1206 14 : AlwaysAssert(image, AipsError);
1207 14 : logIO() << LogOrigin("GridFT", "getImage") << LogIO::NORMAL;
1208 :
1209 14 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
1210 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
1211 : }
1212 14 : weights.resize(sumWeight.shape());
1213 :
1214 14 : convertArray(weights, sumWeight);
1215 : // If the weights are all zero then we cannot normalize
1216 : // otherwise we don't care.
1217 14 : 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 14 : if(useDoubleGrid_p)
1244 : {
1245 28 : 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 : }
1258 : else{
1259 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1260 : //storeArrayAsImage(String("cgrid_ft.im"), image->coordinates(), griddedData);
1261 0 : lattice=arrayLattice;
1262 0 : LatticeFFT::cfft2d(*lattice,false);
1263 : }
1264 :
1265 :
1266 :
1267 : {
1268 14 : Int inx = lattice->shape()(0);
1269 14 : Int iny = lattice->shape()(1);
1270 28 : Vector<Complex> correction(inx);
1271 14 : correction=Complex(1.0, 0.0);
1272 : // Do the Grid-correction
1273 28 : IPosition cursorShape(4, inx, 1, 1, 1);
1274 28 : IPosition axisPath(4, 0, 1, 2, 3);
1275 28 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1276 28 : LatticeIterator<Complex> lix(*lattice, lsx);
1277 2734 : for(lix.reset();!lix.atEnd();lix++) {
1278 2720 : Int pol=lix.position()(2);
1279 2720 : Int chan=lix.position()(3);
1280 2720 : if(weights(pol, chan)!=0.0) {
1281 2720 : gridder->correctX1D(correction, lix.position()(1));
1282 2720 : lix.rwVectorCursor()/=correction;
1283 2720 : if(normalize) {
1284 0 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
1285 0 : 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 : }
1297 :
1298 14 : if(!isTiled) {
1299 : // Check the section from the image BEFORE converting to a lattice
1300 42 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1301 28 : IPosition stride(4, 1);
1302 42 : IPosition trc(blc+image->shape()-stride);
1303 : // Do the copy
1304 14 : IPosition start(4, 0);
1305 14 : image->put(griddedData(blc, trc));
1306 : }
1307 : }
1308 :
1309 14 : image->flush();
1310 14 : image->clearCache();
1311 :
1312 14 : if(!arrayLattice.null()) arrayLattice=0;
1313 14 : if(!lattice.null()) lattice=0;
1314 14 : griddedData.resize();
1315 :
1316 14 : 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 28 : const IPosition latticeShape = weightImage.shape();
1329 :
1330 14 : Int nx=latticeShape(0);
1331 14 : Int ny=latticeShape(1);
1332 :
1333 28 : IPosition loc(2, 0);
1334 28 : IPosition cursorShape(4, nx, ny, 1, 1);
1335 28 : IPosition axisPath(4, 0, 1, 2, 3);
1336 28 : LatticeStepper lsx(latticeShape, cursorShape, axisPath);
1337 28 : 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 2 : 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 2 : Bool retval = true;
1353 2 : Float elpadd=padding_p;
1354 2 : if(toVis_p && withImage)
1355 2 : elpadd=1.0;
1356 : //save the base class variables
1357 2 : 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 2 : outRec.define("cachesize", Int64(cachesize));
1364 2 : outRec.define("tilesize", tilesize);
1365 2 : outRec.define("convtype", convType);
1366 2 : outRec.define("uvscale", uvScale);
1367 2 : outRec.define("uvoffset", uvOffset);
1368 2 : outRec.define("istiled", isTiled);
1369 2 : outRec.define("maxabsdata", maxAbsData);
1370 4 : Vector<Int> center_loc(4), offset_loc(4);
1371 10 : for (Int k=0; k<4 ; k++){
1372 8 : center_loc(k)=centerLoc(k);
1373 8 : offset_loc(k)=offsetLoc(k);
1374 : }
1375 2 : outRec.define("centerloc", center_loc);
1376 2 : outRec.define("offsetloc", offset_loc);
1377 2 : outRec.define("padding", elpadd);
1378 2 : outRec.define("usezero", usezero_p);
1379 2 : outRec.define("nopadding", noPadding_p);
1380 2 : return retval;
1381 : }
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 : }
1439 :
1440 1326 : void GridFT::ok() {
1441 1326 : AlwaysAssert(image, AipsError);
1442 1326 : }
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 7069 : String GridFT::name() const {
1514 :
1515 7069 : return machineName_p;
1516 :
1517 :
1518 : }
1519 :
1520 : } //# NAMESPACE CASA - END
1521 :
|
