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