Line data Source code
1 : //# GridFT.cc: Implementation of GridFT class
2 : //# Copyright (C) 1997-2014
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <msvis/MSVis/VisibilityIterator.h>
29 : #include <casacore/casa/Quanta/UnitMap.h>
30 : #include <casacore/casa/Quanta/UnitVal.h>
31 : #include <casacore/measures/Measures/Stokes.h>
32 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
33 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
34 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
35 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
36 : #include <casacore/coordinates/Coordinates/Projection.h>
37 : #include <casacore/ms/MeasurementSets/MSColumns.h>
38 : #include <casacore/casa/BasicSL/Constants.h>
39 : #include <casacore/scimath/Mathematics/FFTServer.h>
40 : #include <synthesis/TransformMachines/GridFT.h>
41 : #include <casacore/scimath/Mathematics/RigidVector.h>
42 : #include <msvis/MSVis/StokesVector.h>
43 : #include <synthesis/TransformMachines/StokesImageUtil.h>
44 : #include <synthesis/TransformMachines/Utils.h>
45 : #include <msvis/MSVis/VisBuffer.h>
46 : #include <msvis/MSVis/VisSet.h>
47 : #include <casacore/images/Images/ImageInterface.h>
48 : #include <casacore/images/Images/PagedImage.h>
49 : #include <casacore/casa/Containers/Block.h>
50 : #include <casacore/casa/Containers/Record.h>
51 : #include <casacore/casa/Arrays/ArrayLogical.h>
52 : #include <casacore/casa/Arrays/ArrayMath.h>
53 : #include <casacore/casa/Arrays/Array.h>
54 : #include <casacore/casa/Arrays/MaskedArray.h>
55 : #include <casacore/casa/Arrays/Vector.h>
56 : #include <casacore/casa/Arrays/Matrix.h>
57 : #include <casacore/casa/Arrays/Cube.h>
58 : #include <casacore/casa/Arrays/MatrixIter.h>
59 : #include <casacore/casa/BasicSL/String.h>
60 : #include <casacore/casa/Utilities/Assert.h>
61 : #include <casacore/casa/Exceptions/Error.h>
62 : #include <casacore/lattices/Lattices/ArrayLattice.h>
63 : #include <casacore/measures/Measures/UVWMachine.h>
64 : #include <casacore/lattices/Lattices/SubLattice.h>
65 : #include <casacore/lattices/LRegions/LCBox.h>
66 : #include <casacore/lattices/Lattices/LatticeCache.h>
67 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
68 : #include <casacore/lattices/Lattices/LatticeIterator.h>
69 : #include <casacore/lattices/Lattices/LatticeStepper.h>
70 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
71 : #include <casacore/casa/Utilities/CompositeNumber.h>
72 : #include <casacore/casa/OS/Timer.h>
73 : #include <sstream>
74 : #ifdef _OPENMP
75 : #include <omp.h>
76 : #endif
77 :
78 : using namespace casacore;
79 : namespace casa { //# NAMESPACE CASA - BEGIN
80 : using namespace casacore;
81 : // using namespace casa::async;
82 8 : GridFT::GridFT() : FTMachine(), padding_p(1.0), imageCache(0), cachesize(1000000), tilesize(1000), gridder(0), isTiled(false), convType("SF"),
83 8 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
84 8 : usezero_p(false), noPadding_p(false), usePut2_p(false),
85 16 : 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 0 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize), tilesize(itilesize),
91 0 : gridder(0), isTiled(false), convType(iconvType),
92 0 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
93 0 : usezero_p(usezero), noPadding_p(false), usePut2_p(false),
94 0 : machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
95 : {
96 0 : useDoubleGrid_p=useDoublePrec;
97 : // peek=NULL;
98 0 : }
99 :
100 144 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
101 144 : MPosition mLocation, Float padding, Bool usezero, Bool useDoublePrec)
102 144 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
103 144 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
104 144 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
105 288 : 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 144 : mLocation_p=mLocation;
108 144 : tangentSpecified_p=false;
109 144 : useDoubleGrid_p=useDoublePrec;
110 : // peek=NULL;
111 144 : }
112 :
113 0 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
114 0 : MDirection mTangent, Float padding, Bool usezero, Bool useDoublePrec)
115 0 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
116 0 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
117 0 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
118 0 : usePut2_p(false), machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
119 : {
120 0 : mTangent_p=mTangent;
121 0 : tangentSpecified_p=true;
122 0 : useDoubleGrid_p=useDoublePrec;
123 : // peek=NULL;
124 0 : }
125 :
126 4 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
127 : MPosition mLocation, MDirection mTangent, Float padding,
128 4 : Bool usezero, Bool useDoublePrec)
129 4 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
130 4 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
131 4 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
132 8 : 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 70 : void GridFT::init() {
215 :
216 70 : logIO() << LogOrigin("GridFT", "init") << LogIO::NORMAL;
217 :
218 70 : 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 70 : isTiled=false;
239 70 : if(!noPadding_p){
240 70 : CompositeNumber cn(uInt(image->shape()(0)*2));
241 70 : nx = cn.nextLargerEven(Int(padding_p*Float(image->shape()(0))-0.5));
242 70 : ny = cn.nextLargerEven(Int(padding_p*Float(image->shape()(1))-0.5));
243 70 : }
244 : else{
245 0 : nx = image->shape()(0);
246 0 : ny = image->shape()(1);
247 : }
248 70 : npol = image->shape()(2);
249 70 : nchan = image->shape()(3);
250 : // }
251 :
252 70 : sumWeight.resize(npol, nchan);
253 :
254 70 : uvScale.resize(2);
255 70 : uvScale(0)=(Float(nx)*image->coordinates().increment()(0));
256 70 : uvScale(1)=(Float(ny)*image->coordinates().increment()(1));
257 70 : uvOffset.resize(2);
258 70 : uvOffset(0)=nx/2;
259 70 : uvOffset(1)=ny/2;
260 :
261 : // Now set up the gridder. The possibilities are BOX and SF
262 70 : if(gridder) delete gridder; gridder=0;
263 70 : convType.upcase();
264 140 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
265 70 : uvScale, uvOffset,
266 70 : convType);
267 :
268 :
269 70 : convFunc_p.resize();
270 70 : convSupport_p=gridder->cSupport()(0);
271 70 : convSampling_p=gridder->cSampling();
272 70 : 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 70 : if(imageCache) delete imageCache; imageCache=0;
278 :
279 70 : if(isTiled) {
280 0 : Float tileOverlap=0.5;
281 0 : if(convType=="BOX") {
282 0 : tileOverlap=0.0;
283 : }
284 : else {
285 0 : tileOverlap=0.5;
286 0 : tilesize=max(12,tilesize);
287 : }
288 0 : IPosition tileShape=IPosition(4,tilesize,tilesize,npol,nchan);
289 0 : Vector<Float> tileOverlapVec(4);
290 0 : tileOverlapVec=0.0;
291 0 : tileOverlapVec(0)=tileOverlap;
292 0 : tileOverlapVec(1)=tileOverlap;
293 0 : Int tmpCacheVal=static_cast<Int>(cachesize);
294 0 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
295 : tileOverlapVec,
296 0 : (tileOverlap>0.0));
297 :
298 0 : }
299 70 : }
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 55 : 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 55 : IPosition stride(4, 1);
359 110 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
360 110 : IPosition trc(blc+image->shape()-stride);
361 :
362 55 : 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 55 : Vector<Complex> correction(nx);
374 55 : correction=Complex(1.0, 0.0);
375 : // Do the Grid-correction
376 55 : IPosition cursorShape(4, nx, 1, 1, 1);
377 55 : IPosition axisPath(4, 0, 1, 2, 3);
378 55 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
379 55 : 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 55 : }
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 0 : }
416 12 : }
417 :
418 :
419 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
420 : // grid.
421 15 : void GridFT::initializeToSky(ImageInterface<Complex>& iimage,
422 : Matrix<Float>& weight, const VisBuffer& vb)
423 : {
424 : // image always points to the image
425 15 : image=&iimage;
426 15 : toVis_p=false;
427 :
428 15 : init();
429 :
430 : // Initialize the maps for polarization and channel. These maps
431 : // translate visibility indices into image indices
432 15 : initMaps(vb);
433 :
434 :
435 :
436 15 : sumWeight=0.0;
437 15 : weight.resize(sumWeight.shape());
438 15 : weight=0.0;
439 :
440 : // Initialize for in memory or to disk gridding. lattice will
441 : // point to the appropriate Lattice, either the ArrayLattice for
442 : // in memory gridding or to the image for to disk gridding.
443 15 : if(isTiled) {
444 0 : imageCache->flush();
445 0 : image->set(Complex(0.0));
446 0 : lattice=CountedPtr<Lattice<Complex> >(image, false);
447 : }
448 : else {
449 15 : IPosition gridShape(4, nx, ny, npol, nchan);
450 15 : if( !useDoubleGrid_p )
451 : {
452 1 : griddedData.resize(gridShape);
453 1 : griddedData=Complex(0.0);
454 : }
455 15 : if(useDoubleGrid_p){
456 : //griddedData.resize();
457 14 : griddedData2.resize(gridShape);
458 14 : griddedData2=DComplex(0.0);
459 : }
460 15 : image->clearCache();
461 : //iimage.get(griddedData, false);
462 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
463 :
464 15 : }
465 : //AlwaysAssert(lattice, AipsError);
466 15 : }
467 :
468 :
469 :
470 15 : void GridFT::finalizeToSky()
471 : {
472 : //AlwaysAssert(lattice, AipsError);
473 : // Now we flush the cache and report statistics
474 : // For memory based, we don't write anything out yet.
475 : //cerr <<"Time to massage data " << timemass_p << endl;
476 : //cerr <<"Time to grid data " << timegrid_p << endl;
477 15 : timemass_p=0.0;
478 15 : timegrid_p=0.0;
479 15 : if(isTiled) {
480 0 : logIO() << LogOrigin("GridFT", "finalizeToSky") << LogIO::NORMAL;
481 :
482 0 : AlwaysAssert(image, AipsError);
483 0 : AlwaysAssert(imageCache, AipsError);
484 0 : imageCache->flush();
485 0 : ostringstream o;
486 0 : imageCache->showCacheStatistics(o);
487 0 : logIO() << o.str() << LogIO::POST;
488 0 : }
489 : // peek->terminate();
490 : // peek->reset();
491 15 : }
492 :
493 :
494 :
495 0 : Array<Complex>* GridFT::getDataPointer(const IPosition& centerLoc2D,
496 : Bool readonly) {
497 : Array<Complex>* result;
498 : // Is tiled: get tiles and set up offsets
499 0 : centerLoc(0)=centerLoc2D(0);
500 0 : centerLoc(1)=centerLoc2D(1);
501 0 : result=&imageCache->tile(offsetLoc,centerLoc, readonly);
502 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
503 0 : return result;
504 : }
505 :
506 : #define NEED_UNDERSCORES
507 : #if defined(NEED_UNDERSCORES)
508 : #define ggrid ggrid_
509 : #define dgrid dgrid_
510 : #define ggrids ggrids_
511 : #define sectggridd sectggridd_
512 : #define sectggrids sectggrids_
513 : #define sectdgrid sectdgrid_
514 : #define locuvw locuvw_
515 : #endif
516 :
517 : extern "C" {
518 : void locuvw(const Double*, const Double*, const Double*, const Int*, const Double*, const Double*, const Int*,
519 : Int*, Int*, Complex*, const Int*, const Int*, const Double*);
520 : void ggrid(Double*,
521 : Double*,
522 : const Complex*,
523 : Int*,
524 : Int*,
525 : Int*,
526 : const Int*,
527 : const Int*,
528 : const Float*,
529 : Int*,
530 : Int*,
531 : Double*,
532 : Double*,
533 : DComplex*,
534 : Int*,
535 : Int*,
536 : Int *,
537 : Int *,
538 : const Double*,
539 : const Double*,
540 : Int*,
541 : Int*,
542 : Double*,
543 : Int*,
544 : Int*,
545 : Double*);
546 : void sectggridd(const Complex*,
547 : const Int*,
548 : const Int*,
549 : const Int*,
550 : const Int*,
551 : const Int*,
552 : const Float*,
553 : const Int*,
554 : DComplex*,
555 : const Int*,
556 : const Int*,
557 : const Int *,
558 : const Int *,
559 : //support
560 : const Int*,
561 : const Int*,
562 : const Double*,
563 : const Int*,
564 : const Int*,
565 : Double*,
566 : const Int*,
567 : const Int*,
568 : const Int*,
569 : const Int*,
570 : const Int*,
571 : const Int*,
572 : //loc, off, phasor
573 : const Int*,
574 : const Int*,
575 : const Complex*);
576 : //single precision gridder
577 : void sectggrids(const Complex*,
578 : const Int*,
579 : const Int*,
580 : const Int*,
581 : const Int*,
582 : const Int*,
583 : const Float*,
584 : const Int*,
585 : Complex*,
586 : const Int*,
587 : const Int*,
588 : const Int *,
589 : const Int *,
590 : const Int*,
591 : const Int*,
592 : const Double*,
593 : const Int*,
594 : const Int*,
595 : Double*,
596 : const Int*,
597 : const Int*,
598 : const Int*,
599 : const Int*,
600 : const Int*,
601 : const Int*,
602 : //loc, off, phasor
603 : const Int*,
604 : const Int*,
605 : const Complex*);
606 : void ggrids(Double*,
607 : Double*,
608 : const Complex*,
609 : Int*,
610 : Int*,
611 : Int*,
612 : const Int*,
613 : const Int*,
614 : const Float*,
615 : Int*,
616 : Int*,
617 : Double*,
618 : Double*,
619 : Complex*,
620 : Int*,
621 : Int*,
622 : Int *,
623 : Int *,
624 : const Double*,
625 : const Double*,
626 : Int*,
627 : Int*,
628 : Double*,
629 : Int*,
630 : Int*,
631 : Double*);
632 :
633 : void dgrid(Double*,
634 : Double*,
635 : Complex*,
636 : Int*,
637 : Int*,
638 : const Int*,
639 : const Int*,
640 : Int*,
641 : Int*,
642 : Double*,
643 : Double*,
644 : const Complex*,
645 : Int*,
646 : Int*,
647 : Int *,
648 : Int *,
649 : const Double*,
650 : const Double*,
651 : Int*,
652 : Int*,
653 : Double*,
654 : Int*,
655 : Int*);
656 :
657 : void sectdgrid(Complex*,
658 : const Int*,
659 : const Int*,
660 : const Int*,
661 : const Int*,
662 : const Int*,
663 : const Complex*,
664 : const Int*,
665 : const Int*,
666 : const Int *,
667 : const Int *,
668 : //support
669 : const Int*,
670 : const Int*,
671 : const Double*,
672 : const Int*,
673 : const Int*,
674 : //rbeg, rend, loc, off, phasor
675 : const Int*,
676 : const Int*,
677 : const Int*,
678 : const Int*,
679 : const Complex*);
680 : }
681 463 : void GridFT::put(const VisBuffer& vb, Int row, Bool dopsf,
682 : FTMachine::Type type)
683 : {
684 :
685 463 : gridOk(convSupport_p);
686 : // peek->setVBPtr(&vb);
687 :
688 : //Check if ms has changed then cache new spw and chan selection
689 463 : if(vb.newMS())
690 0 : matchAllSpwChans(vb);
691 :
692 : //Here we redo the match or use previous match
693 :
694 : //Channel matching for the actual spectral window of buffer
695 463 : if(doConversion_p[vb.spectralWindow()]){
696 0 : matchChannel(vb.spectralWindow(), vb);
697 : }
698 : else{
699 463 : chanMap.resize();
700 463 : chanMap=multiChanMap_p[vb.spectralWindow()];
701 : }
702 :
703 : //No point in reading data if its not matching in frequency
704 463 : if(max(chanMap)==-1)
705 : {
706 : // peek->reset();
707 0 : return;
708 : }
709 :
710 463 : Timer tim;
711 463 : tim.mark();
712 :
713 : const Matrix<Float> *imagingweight;
714 463 : imagingweight=&(vb.imagingWeight());
715 :
716 463 : if(dopsf) {type=FTMachine::PSF;}
717 :
718 463 : Cube<Complex> data;
719 : //Fortran gridder need the flag as ints
720 463 : Cube<Int> flags;
721 463 : Matrix<Float> elWeight;
722 463 : interpolateFrequencyTogrid(vb, *imagingweight,data, flags, elWeight, type);
723 :
724 :
725 : Bool iswgtCopy;
726 : const Float *wgtStorage;
727 463 : wgtStorage=elWeight.getStorage(iswgtCopy);
728 :
729 : Bool isCopy;
730 : const Complex *datStorage;
731 463 : if(!dopsf)
732 200 : datStorage=data.getStorage(isCopy);
733 : else
734 263 : datStorage=0;
735 : // If row is -1 then we pass through all rows
736 : Int startRow, endRow, nRow;
737 463 : if (row==-1) {
738 463 : nRow=vb.nRow();
739 463 : startRow=0;
740 463 : endRow=nRow-1;
741 : } else {
742 0 : nRow=1;
743 0 : startRow=row;
744 0 : endRow=row;
745 : }
746 :
747 : //cerr << "nRow " << nRow << endl;
748 :
749 :
750 : // Get the uvws in a form that Fortran can use and do that
751 : // necessary phase rotation. On a Pentium Pro 200 MHz
752 : // when null, this step takes about 50us per uvw point. This
753 : // is just barely noticeable for Stokes I continuum and
754 : // irrelevant for other cases.
755 :
756 : Bool del;
757 : // IPosition s(flags.shape());
758 463 : const IPosition &fs=flags.shape();
759 463 : std::vector<Int> s(fs.begin(),fs.end());
760 463 : Int nvispol=s[0];
761 463 : Int nvischan=s[1];
762 463 : Int nvisrow=s[2];
763 463 : Matrix<Double> uvw(3, vb.uvw().nelements());
764 463 : uvw=0.0;
765 463 : Vector<Double> dphase(vb.uvw().nelements());
766 463 : Cube<Int> loc(2, nvischan, nRow);
767 463 : Cube<Int> off(2, nvischan, nRow);
768 463 : Matrix<Complex> phasor(nvischan, nRow);
769 463 : Int csamp=convSampling_p;
770 463 : dphase=0.0;
771 : //NEGATING to correct for an image inversion problem
772 320178 : for (Int i=startRow;i<=endRow;i++) {
773 959145 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
774 319715 : uvw(2,i)=vb.uvw()(i)(2);
775 : }
776 463 : timemass_p +=tim.real();
777 463 : tim.mark();
778 463 : rotateUVW(uvw, dphase, vb);
779 463 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
780 : Bool delphase;
781 463 : Complex * phasorstor=phasor.getStorage(delphase);
782 463 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
783 463 : const Double * scalestor=uvScale.getStorage(del);
784 463 : const Double * offsetstor=uvOffset.getStorage(del);
785 463 : const Double* uvstor= uvw.getStorage(del);
786 463 : Int * locstor=loc.getStorage(del);
787 463 : Int * offstor=off.getStorage(del);
788 463 : const Double *dpstor=dphase.getStorage(del);
789 : //Vector<Double> f1=interpVisFreq_p.copy();
790 463 : Int nvchan=nvischan;
791 : Int irow;
792 463 : Double cinv=Double(1.0)/C::c;
793 463 : Int dow=0;
794 463 : Int nth=1;
795 : #ifdef _OPENMP
796 463 : if(numthreads_p >0){
797 0 : nth=min(numthreads_p, omp_get_max_threads());
798 : }
799 : else{
800 463 : nth= omp_get_max_threads();
801 : }
802 463 : nth=min(4,nth);
803 : #endif
804 :
805 :
806 463 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvchan, scalestor, offsetstor, csamp, phasorstor, uvstor, locstor, offstor, dpstor, cinv, dow) shared(startRow, endRow) num_threads(nth)
807 : {
808 : #pragma omp for
809 : for (irow=startRow; irow<=endRow; ++irow){
810 : //locateuvw(uvstor,dpstor, visfreqstor, nvchan, scalestor, offsetstor, csamp,
811 : // locstor,
812 : // offstor, phasorstor, irow);
813 : locuvw(uvstor, dpstor, visfreqstor, &nvchan, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
814 : }
815 :
816 : }//end pragma parallel
817 : // Take care of translation of Bools to Integer
818 : //////TESTOO
819 : //ofstream myfile;
820 : //myfile.open ("putLoc.txt", ios::out | ios::app | ios::ate );
821 : //myfile << vb.rowIds()(0) << " uv " << uvw.column(0) << " loc " << loc(0,0,0) << ", " << loc(1,0,0) << "\n" << endl;
822 : //myfile.close();
823 : ///////////////
824 :
825 463 : Int idopsf=0;
826 463 : if(dopsf) idopsf=1;
827 :
828 :
829 463 : Vector<Int> rowFlags(vb.nRow());
830 463 : rowFlags=0;
831 463 : rowFlags(vb.flagRow())=true;
832 463 : if(!usezero_p) {
833 320178 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
834 319715 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
835 : }
836 : }
837 :
838 : /////////////Some extra stuff for openmp
839 :
840 : Int x0, y0, nxsub, nysub, ixsub, iysub, icounter, ix, iy;
841 463 : Int csupp=convSupport_p;
842 :
843 463 : const Double * convfuncstor=(gridder->cFunction()).getStorage(del);
844 :
845 : // cerr <<"Poffset " << min(off) << " " << max(off) << " length " << gridder->cFunction().shape() << endl;
846 :
847 :
848 :
849 463 : ixsub=1;
850 463 : iysub=1;
851 463 : if (nth >3){
852 463 : ixsub=2;
853 463 : iysub=2;
854 : }
855 0 : else if(nth >1){
856 0 : ixsub=2;
857 0 : iysub=1;
858 : }
859 :
860 463 : x0=1;
861 463 : y0=1;
862 463 : nxsub=nx;
863 463 : nysub=ny;
864 463 : Int rbeg=startRow+1;
865 463 : Int rend=endRow+1;
866 463 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
867 2315 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
868 1852 : sumwgt[icounter].resize(sumWeight.shape());
869 1852 : sumwgt[icounter].set(0.0);
870 : }
871 463 : const Int* pmapstor=polMap.getStorage(del);
872 463 : const Int *cmapstor=chanMap.getStorage(del);
873 463 : Int nc=nchan;
874 463 : Int np=npol;
875 463 : Int nxp=nx;
876 463 : Int nyp=ny;
877 463 : const Int * flagstor=flags.getStorage(del);
878 463 : const Int * rowflagstor=rowFlags.getStorage(del);
879 : ////////////////////////
880 :
881 : Bool gridcopy;
882 463 : if(useDoubleGrid_p){
883 263 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
884 263 : #pragma omp parallel default(none) private(icounter,ix,iy,x0,y0,nxsub,nysub, del) firstprivate(idopsf, datStorage, wgtStorage, flagstor, rowflagstor, convfuncstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc,ixsub, iysub, rend, rbeg, csamp, csupp, nvispol, nvischan, nvisrow, phasorstor, locstor, offstor) shared(sumwgt) num_threads(ixsub*iysub)
885 :
886 : {
887 : //cerr << "numthreads " << omp_get_num_threads() << endl;
888 : #pragma omp for
889 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
890 : //cerr << "thread id " << omp_get_thread_num() << endl;
891 : ix= (icounter+1)-((icounter)/ixsub)*ixsub;
892 : iy=(icounter)/ixsub+1;
893 : y0=(nyp/iysub)*(iy-1)+1;
894 : nysub=nyp/iysub;
895 : if( iy == iysub) {
896 : nysub=nyp-(nyp/iysub)*(iy-1);
897 : }
898 : x0=(nxp/ixsub)*(ix-1)+1;
899 : nxsub=nxp/ixsub;
900 : if( ix == ixsub){
901 : nxsub=nxp-(nxp/ixsub)*(ix-1);
902 : }
903 : sectggridd(datStorage,
904 : &nvispol,
905 : &nvischan,
906 : &idopsf,
907 : flagstor,
908 : rowflagstor,
909 : wgtStorage,
910 : &nvisrow,
911 : gridstor,
912 : &nxp,
913 : &nyp,
914 : &np,
915 : &nc,
916 : &csupp,
917 : &csamp,
918 : convfuncstor,
919 : cmapstor,
920 : pmapstor,
921 : (sumwgt[icounter]).getStorage(del),
922 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
923 : phasorstor);
924 : }//end for
925 : }// end pragma parallel
926 1315 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
927 1052 : sumWeight=sumWeight+sumwgt[icounter];
928 : }
929 : //phasor.putStorage(phasorstor, delphase);
930 263 : griddedData2.putStorage(gridstor, gridcopy);
931 : }
932 : else{
933 200 : Complex *gridstor=griddedData.getStorage(gridcopy);
934 200 : #pragma omp parallel default(none) private(icounter,ix,iy,x0,y0,nxsub,nysub, del) firstprivate(idopsf, datStorage, wgtStorage, flagstor, rowflagstor, convfuncstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc,ixsub, iysub, rend, rbeg, csamp, csupp, nvispol, nvischan, nvisrow, phasorstor, locstor, offstor) shared(sumwgt) num_threads(ixsub*iysub)
935 : {
936 : //cerr << "numthreads " << omp_get_num_threads() << endl;
937 : #pragma omp for
938 :
939 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
940 : //cerr << "thread id " << omp_get_thread_num() << endl;
941 : ix= (icounter+1)-((icounter)/ixsub)*ixsub;
942 : iy=(icounter)/ixsub+1;
943 : y0=(nyp/iysub)*(iy-1)+1;
944 : nysub=nyp/iysub;
945 : if( iy == iysub) {
946 : nysub=nyp-(nyp/iysub)*(iy-1);
947 : }
948 : x0=(nxp/ixsub)*(ix-1)+1;
949 : nxsub=nxp/ixsub;
950 : if( ix == ixsub){
951 : nxsub=nxp-(nxp/ixsub)*(ix-1);
952 : }
953 : //cerr << "x0 " << x0 << " y0 " << y0 << " nxsub " << nxsub << " nysub " << nysub << endl;
954 : sectggrids(datStorage,
955 : &nvispol,
956 : &nvischan,
957 : &idopsf,
958 : flagstor,
959 : rowflagstor,
960 : wgtStorage,
961 : &nvisrow,
962 : gridstor,
963 : &nxp,
964 : &nyp,
965 : &np,
966 : &nc,
967 : &csupp,
968 : &csamp,
969 : convfuncstor,
970 : cmapstor,
971 : pmapstor,
972 : (sumwgt[icounter]).getStorage(del),
973 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
974 : phasorstor);
975 : }//end for
976 : }// end pragma parallel
977 1000 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
978 800 : sumWeight=sumWeight+sumwgt[icounter];
979 : }
980 :
981 200 : griddedData.putStorage(gridstor, gridcopy);
982 : }
983 : // cerr << "sunweight " << sumWeight << endl;
984 :
985 463 : timegrid_p+=tim.real();
986 :
987 463 : if(!dopsf)
988 200 : data.freeStorage(datStorage, isCopy);
989 463 : elWeight.freeStorage(wgtStorage,iswgtCopy);
990 :
991 : // peek->reset();
992 463 : }
993 :
994 0 : void GridFT::modifyConvFunc(const Vector<Double>& convFunc, Int convSupport, Int convSampling){
995 0 : convFunc_p.resize();
996 0 : convFunc_p=convFunc;
997 0 : convSupport_p=convSupport;
998 0 : convSampling_p=convSampling;
999 :
1000 0 : }
1001 :
1002 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 266 : Cube<Complex> data;
1059 266 : 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 266 : Cube<Int> loc(2, nvc, nvisrow);
1077 266 : Cube<Int> off(2, nvc, nRow);
1078 266 : 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 266 : 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 266 : npart=4;
1157 : }
1158 0 : 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 266 : }
1197 :
1198 :
1199 :
1200 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1201 : // return the resulting image
1202 15 : ImageInterface<Complex>& GridFT::getImage(Matrix<Float>& weights, Bool normalize)
1203 : {
1204 : //AlwaysAssert(lattice, AipsError);
1205 15 : AlwaysAssert(gridder, AipsError);
1206 15 : AlwaysAssert(image, AipsError);
1207 15 : logIO() << LogOrigin("GridFT", "getImage") << LogIO::NORMAL;
1208 :
1209 15 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
1210 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
1211 : }
1212 15 : weights.resize(sumWeight.shape());
1213 :
1214 15 : convertArray(weights, sumWeight);
1215 : // If the weights are all zero then we cannot normalize
1216 : // otherwise we don't care.
1217 15 : if(normalize&&max(weights)==0.0) {
1218 : logIO() << LogIO::SEVERE << "No useful data in GridFT: weights all zero"
1219 0 : << LogIO::POST;
1220 : }
1221 : else {
1222 :
1223 : //const IPosition latticeShape = lattice->shape();
1224 :
1225 : //logIO() << LogIO::DEBUGGING
1226 : // << "Starting FFT and scaling of image" << LogIO::POST;
1227 :
1228 :
1229 :
1230 : // if(useDoubleGrid_p){
1231 : // convertArray(griddedData, griddedData2);
1232 : // //Don't need the double-prec grid anymore...
1233 : // griddedData2.resize();
1234 : // }
1235 :
1236 : // x and y transforms
1237 : // LatticeFFT::cfft2d(*lattice,false);
1238 : //
1239 : // Retain the double precision grid for FFT as well. Convert it
1240 : // to single precision just after (since images are still single
1241 : // precision).
1242 : //
1243 15 : if(useDoubleGrid_p)
1244 : {
1245 14 : ArrayLattice<DComplex> darrayLattice(griddedData2);
1246 :
1247 14 : LatticeFFT::cfft2d(darrayLattice,false);
1248 14 : griddedData.resize(griddedData2.shape());
1249 14 : convertArray(griddedData, griddedData2);
1250 :
1251 14 : SynthesisUtilMethods::getResource("mem peak in getImage");
1252 :
1253 : //Don't need the double-prec grid anymore...
1254 14 : griddedData2.resize();
1255 14 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1256 14 : lattice=arrayLattice;
1257 14 : }
1258 : else{
1259 1 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1260 : //storeArrayAsImage(String("cgrid_ft.im"), image->coordinates(), griddedData);
1261 1 : lattice=arrayLattice;
1262 1 : LatticeFFT::cfft2d(*lattice,false);
1263 : }
1264 :
1265 :
1266 :
1267 : {
1268 15 : Int inx = lattice->shape()(0);
1269 15 : Int iny = lattice->shape()(1);
1270 15 : Vector<Complex> correction(inx);
1271 15 : correction=Complex(1.0, 0.0);
1272 : // Do the Grid-correction
1273 15 : IPosition cursorShape(4, inx, 1, 1, 1);
1274 15 : IPosition axisPath(4, 0, 1, 2, 3);
1275 15 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1276 15 : LatticeIterator<Complex> lix(*lattice, lsx);
1277 2835 : for(lix.reset();!lix.atEnd();lix++) {
1278 2820 : Int pol=lix.position()(2);
1279 2820 : Int chan=lix.position()(3);
1280 2820 : if(weights(pol, chan)!=0.0) {
1281 2820 : gridder->correctX1D(correction, lix.position()(1));
1282 2820 : lix.rwVectorCursor()/=correction;
1283 2820 : if(normalize) {
1284 100 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
1285 100 : lix.rwCursor()*=rnorm;
1286 : }
1287 : else {
1288 2720 : Complex rnorm(Float(inx)*Float(iny));
1289 2720 : lix.rwCursor()*=rnorm;
1290 : }
1291 : }
1292 : else {
1293 0 : lix.woCursor()=0.0;
1294 : }
1295 : }
1296 15 : }
1297 :
1298 15 : if(!isTiled) {
1299 : // Check the section from the image BEFORE converting to a lattice
1300 30 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1301 15 : IPosition stride(4, 1);
1302 30 : IPosition trc(blc+image->shape()-stride);
1303 : // Do the copy
1304 15 : IPosition start(4, 0);
1305 15 : image->put(griddedData(blc, trc));
1306 15 : }
1307 : }
1308 :
1309 15 : image->flush();
1310 15 : image->clearCache();
1311 :
1312 15 : if(!arrayLattice.null()) arrayLattice=0;
1313 15 : if(!lattice.null()) lattice=0;
1314 15 : griddedData.resize();
1315 :
1316 15 : return *image;
1317 : }
1318 :
1319 : // Get weight image
1320 14 : void GridFT::getWeightImage(ImageInterface<Float>& weightImage, Matrix<Float>& weights)
1321 : {
1322 :
1323 14 : logIO() << LogOrigin("GridFT", "getWeightImage") << LogIO::NORMAL;
1324 :
1325 14 : weights.resize(sumWeight.shape());
1326 14 : convertArray(weights,sumWeight);
1327 :
1328 14 : const IPosition latticeShape = weightImage.shape();
1329 :
1330 14 : Int nx=latticeShape(0);
1331 14 : Int ny=latticeShape(1);
1332 :
1333 14 : IPosition loc(2, 0);
1334 14 : IPosition cursorShape(4, nx, ny, 1, 1);
1335 14 : IPosition axisPath(4, 0, 1, 2, 3);
1336 14 : LatticeStepper lsx(latticeShape, cursorShape, axisPath);
1337 14 : LatticeIterator<Float> lix(weightImage, lsx);
1338 28 : for(lix.reset();!lix.atEnd();lix++) {
1339 14 : Int pol=lix.position()(2);
1340 14 : Int chan=lix.position()(3);
1341 14 : lix.rwCursor()=weights(pol,chan);
1342 : }
1343 14 : }
1344 :
1345 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 2 : 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 2 : }
1382 :
1383 0 : Bool GridFT::fromRecord(String& error,
1384 : const RecordInterface& inRec)
1385 : {
1386 0 : Bool retval = true;
1387 0 : if(!FTMachine::fromRecord(error, inRec))
1388 0 : return false;
1389 0 : gridder=0; imageCache=0; lattice=0; arrayLattice=0;
1390 : //For some reason int64 seems to be getting lost...
1391 : //this is not an important parameter...so filing a jira
1392 0 : if(inRec.isDefined("cachesize"))
1393 0 : cachesize=inRec.asInt64("cachesize");
1394 : else
1395 0 : cachesize=1000000;
1396 0 : inRec.get("tilesize", tilesize);
1397 0 : inRec.get("convtype", convType);
1398 0 : inRec.get("uvscale", uvScale);
1399 0 : inRec.get("uvoffset", uvOffset);
1400 0 : inRec.get("istiled", isTiled);
1401 0 : inRec.get("maxabsdata", maxAbsData);
1402 0 : Vector<Int> center_loc(4), offset_loc(4);
1403 0 : inRec.get("centerloc", center_loc);
1404 0 : inRec.get("offsetloc", offset_loc);
1405 0 : uInt ndim4 = 4;
1406 0 : centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2),
1407 0 : center_loc(3));
1408 0 : offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2),
1409 0 : offset_loc(3));
1410 0 : inRec.get("padding", padding_p);
1411 0 : inRec.get("usezero", usezero_p);
1412 0 : inRec.get("nopadding", noPadding_p);
1413 :
1414 0 : machineName_p="GridFT";
1415 : ///setup some of the parameters
1416 0 : init();
1417 : /*if(!cmplxImage_p.null()){
1418 : //FTMachine::fromRecord would have recovered the image
1419 : // Might be changing the shape of sumWeight
1420 :
1421 : if(isTiled) {
1422 : lattice=CountedPtr<Lattice<Complex> >(image, false);
1423 : }
1424 : else {
1425 : // Make the grid the correct shape and turn it into an array lattice
1426 : // Check the section from the image BEFORE converting to a lattice
1427 : if(!toVis_p){
1428 : IPosition gridShape(4, nx, ny, npol, nchan);
1429 : griddedData.resize(gridShape);
1430 : griddedData=Complex(0.0);
1431 : }
1432 : else{
1433 : prepGridForDegrid();
1434 : }
1435 : }
1436 : };*/
1437 0 : return retval;
1438 0 : }
1439 :
1440 1727 : void GridFT::ok() {
1441 1727 : AlwaysAssert(image, AipsError);
1442 1727 : }
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 :
|
