Line data Source code
1 : //# GridFT.cc: Implementation of GridFT class
2 : //# Copyright (C) 1997-2014
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <msvis/MSVis/VisibilityIterator.h>
29 : #include <casacore/casa/Quanta/UnitMap.h>
30 : #include <casacore/casa/Quanta/UnitVal.h>
31 : #include <casacore/measures/Measures/Stokes.h>
32 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
33 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
34 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
35 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
36 : #include <casacore/coordinates/Coordinates/Projection.h>
37 : #include <casacore/ms/MeasurementSets/MSColumns.h>
38 : #include <casacore/casa/BasicSL/Constants.h>
39 : #include <casacore/scimath/Mathematics/FFTServer.h>
40 : #include <synthesis/TransformMachines/GridFT.h>
41 : #include <casacore/scimath/Mathematics/RigidVector.h>
42 : #include <msvis/MSVis/StokesVector.h>
43 : #include <synthesis/TransformMachines/StokesImageUtil.h>
44 : #include <synthesis/TransformMachines/Utils.h>
45 : #include <msvis/MSVis/VisBuffer.h>
46 : #include <msvis/MSVis/VisSet.h>
47 : #include <casacore/images/Images/ImageInterface.h>
48 : #include <casacore/images/Images/PagedImage.h>
49 : #include <casacore/casa/Containers/Block.h>
50 : #include <casacore/casa/Containers/Record.h>
51 : #include <casacore/casa/Arrays/ArrayLogical.h>
52 : #include <casacore/casa/Arrays/ArrayMath.h>
53 : #include <casacore/casa/Arrays/Array.h>
54 : #include <casacore/casa/Arrays/MaskedArray.h>
55 : #include <casacore/casa/Arrays/Vector.h>
56 : #include <casacore/casa/Arrays/Matrix.h>
57 : #include <casacore/casa/Arrays/Cube.h>
58 : #include <casacore/casa/Arrays/MatrixIter.h>
59 : #include <casacore/casa/BasicSL/String.h>
60 : #include <casacore/casa/Utilities/Assert.h>
61 : #include <casacore/casa/Exceptions/Error.h>
62 : #include <casacore/lattices/Lattices/ArrayLattice.h>
63 : #include <casacore/measures/Measures/UVWMachine.h>
64 : #include <casacore/lattices/Lattices/SubLattice.h>
65 : #include <casacore/lattices/LRegions/LCBox.h>
66 : #include <casacore/lattices/Lattices/LatticeCache.h>
67 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
68 : #include <casacore/lattices/Lattices/LatticeIterator.h>
69 : #include <casacore/lattices/Lattices/LatticeStepper.h>
70 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
71 : #include <casacore/casa/Utilities/CompositeNumber.h>
72 : #include <casacore/casa/OS/Timer.h>
73 : #include <sstream>
74 : #ifdef _OPENMP
75 : #include <omp.h>
76 : #endif
77 :
78 : using namespace casacore;
79 : namespace casa { //# NAMESPACE CASA - BEGIN
80 : using namespace casacore;
81 : // using namespace casa::async;
82 0 : GridFT::GridFT() : FTMachine(), padding_p(1.0), imageCache(0), cachesize(1000000), tilesize(1000), gridder(0), isTiled(false), convType("SF"),
83 0 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
84 0 : usezero_p(false), noPadding_p(false), usePut2_p(false),
85 0 : machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0) {
86 :
87 0 : }
88 0 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType, Float padding,
89 0 : Bool usezero, Bool useDoublePrec)
90 0 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize), tilesize(itilesize),
91 0 : gridder(0), isTiled(false), convType(iconvType),
92 0 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
93 0 : usezero_p(usezero), noPadding_p(false), usePut2_p(false),
94 0 : machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
95 : {
96 0 : useDoubleGrid_p=useDoublePrec;
97 : // peek=NULL;
98 0 : }
99 :
100 53 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
101 53 : MPosition mLocation, Float padding, Bool usezero, Bool useDoublePrec)
102 53 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
103 53 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
104 53 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
105 106 : 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 53 : mLocation_p=mLocation;
108 53 : tangentSpecified_p=false;
109 53 : useDoubleGrid_p=useDoublePrec;
110 : // peek=NULL;
111 53 : }
112 :
113 0 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
114 0 : MDirection mTangent, Float padding, Bool usezero, Bool useDoublePrec)
115 0 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
116 0 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
117 0 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
118 0 : usePut2_p(false), machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
119 : {
120 0 : mTangent_p=mTangent;
121 0 : tangentSpecified_p=true;
122 0 : useDoubleGrid_p=useDoublePrec;
123 : // peek=NULL;
124 0 : }
125 :
126 0 : GridFT::GridFT(Long icachesize, Int itilesize, String iconvType,
127 : MPosition mLocation, MDirection mTangent, Float padding,
128 0 : Bool usezero, Bool useDoublePrec)
129 0 : : FTMachine(), padding_p(padding), imageCache(0), cachesize(icachesize),
130 0 : tilesize(itilesize), gridder(0), isTiled(false), convType(iconvType), maxAbsData(0.0), centerLoc(IPosition(4,0)),
131 0 : offsetLoc(IPosition(4,0)), usezero_p(usezero), noPadding_p(false),
132 0 : usePut2_p(false),machineName_p("GridFT"), timemass_p(0.0), timegrid_p(0.0), convFunc_p(0), convSampling_p(1), convSupport_p(0)
133 : {
134 0 : mLocation_p=mLocation;
135 0 : mTangent_p=mTangent;
136 0 : tangentSpecified_p=true;
137 0 : useDoubleGrid_p=useDoublePrec;
138 : // peek=NULL;
139 0 : }
140 :
141 0 : GridFT::GridFT(const RecordInterface& stateRec)
142 0 : : FTMachine()
143 : {
144 : // Construct from the input state record
145 0 : String error;
146 0 : if (!fromRecord(error, stateRec)) {
147 0 : throw (AipsError("Failed to create gridder: " + error));
148 : };
149 0 : timemass_p=0.0;
150 0 : timegrid_p=0.0;
151 0 : timedegrid_p=0.0;
152 : // peek=NULL;
153 0 : }
154 :
155 : //----------------------------------------------------------------------
156 106 : GridFT& GridFT::operator=(const GridFT& other)
157 : {
158 106 : if(this!=&other) {
159 : //Do the base parameters
160 106 : FTMachine::operator=(other);
161 :
162 : //private params
163 106 : imageCache=other.imageCache;
164 106 : cachesize=other.cachesize;
165 106 : tilesize=other.tilesize;
166 106 : convType=other.convType;
167 106 : convType.upcase();
168 106 : uvScale.resize();
169 106 : uvOffset.resize();
170 106 : uvScale.assign(other.uvScale);
171 106 : uvOffset.assign(other.uvOffset);
172 106 : if(other.gridder==0)
173 106 : gridder=0;
174 : else{
175 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
176 0 : uvScale, uvOffset,
177 0 : convType);
178 : }
179 106 : convFunc_p.resize();
180 106 : convSupport_p=other.convSupport_p;
181 106 : convSampling_p=other.convSampling_p;
182 106 : convFunc_p=other.convFunc_p;
183 106 : isTiled=other.isTiled;
184 : //lattice=other.lattice;
185 106 : lattice=0;
186 106 : tilesize=other.tilesize;
187 106 : arrayLattice=0;
188 106 : maxAbsData=other.maxAbsData;
189 106 : centerLoc=other.centerLoc;
190 106 : offsetLoc=other.offsetLoc;
191 106 : padding_p=other.padding_p;
192 106 : usezero_p=other.usezero_p;
193 106 : noPadding_p=other.noPadding_p;
194 106 : machineName_p="GridFT";
195 106 : timemass_p=0.0;
196 106 : timegrid_p=0.0;
197 : // peek = other.peek;
198 : };
199 106 : return *this;
200 : };
201 :
202 : //----------------------------------------------------------------------
203 106 : GridFT::GridFT(const GridFT& other) : FTMachine(), machineName_p("GridFT")
204 : {
205 106 : operator=(other);
206 106 : }
207 :
208 : //-----------------------------------------------------------------------
209 0 : FTMachine* GridFT::cloneFTM(){
210 0 : return new GridFT(*this);
211 : }
212 :
213 : //----------------------------------------------------------------------
214 5 : void GridFT::init() {
215 :
216 5 : logIO() << LogOrigin("GridFT", "init") << LogIO::NORMAL;
217 :
218 5 : 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 5 : isTiled=false;
239 5 : if(!noPadding_p){
240 5 : CompositeNumber cn(uInt(image->shape()(0)*2));
241 5 : nx = cn.nextLargerEven(Int(padding_p*Float(image->shape()(0))-0.5));
242 5 : ny = cn.nextLargerEven(Int(padding_p*Float(image->shape()(1))-0.5));
243 5 : }
244 : else{
245 0 : nx = image->shape()(0);
246 0 : ny = image->shape()(1);
247 : }
248 5 : npol = image->shape()(2);
249 5 : nchan = image->shape()(3);
250 : // }
251 :
252 5 : sumWeight.resize(npol, nchan);
253 :
254 5 : uvScale.resize(2);
255 5 : uvScale(0)=(Float(nx)*image->coordinates().increment()(0));
256 5 : uvScale(1)=(Float(ny)*image->coordinates().increment()(1));
257 5 : uvOffset.resize(2);
258 5 : uvOffset(0)=nx/2;
259 5 : uvOffset(1)=ny/2;
260 :
261 : // Now set up the gridder. The possibilities are BOX and SF
262 5 : if(gridder) delete gridder; gridder=0;
263 5 : convType.upcase();
264 10 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
265 5 : uvScale, uvOffset,
266 5 : convType);
267 :
268 :
269 5 : convFunc_p.resize();
270 5 : convSupport_p=gridder->cSupport()(0);
271 5 : convSampling_p=gridder->cSampling();
272 5 : 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 5 : if(imageCache) delete imageCache; imageCache=0;
278 :
279 5 : 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 5 : }
300 :
301 : // This is nasty, we should use CountedPointers here.
302 316 : GridFT::~GridFT() {
303 158 : if(imageCache) delete imageCache; imageCache=0;
304 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
305 158 : if(gridder) delete gridder; gridder=0;
306 316 : }
307 :
308 : // Initialize for a transform from the Sky domain. This means that
309 : // we grid-correct, and FFT the image
310 :
311 4 : void GridFT::initializeToVis(ImageInterface<Complex>& iimage,
312 : const VisBuffer& vb)
313 : {
314 4 : image=&iimage;
315 4 : toVis_p=true;
316 :
317 4 : ok();
318 :
319 4 : init();
320 : // peek->reset();
321 : // Initialize the maps for polarization and channel. These maps
322 : // translate visibility indices into image indices
323 4 : 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 4 : prepGridForDegrid();
342 :
343 : //AlwaysAssert(lattice, AipsError);
344 :
345 :
346 :
347 4 : }
348 :
349 4 : void GridFT::prepGridForDegrid(){
350 4 : IPosition gridShape(4, nx, ny, npol, nchan);
351 4 : 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 4 : griddedData.set(Complex(0.0));
357 :
358 4 : IPosition stride(4, 1);
359 8 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
360 8 : IPosition trc(blc+image->shape()-stride);
361 :
362 4 : IPosition start(4, 0);
363 4 : griddedData(blc, trc) = image->getSlice(start, image->shape());
364 :
365 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
366 4 : arrayLattice = new ArrayLattice<Complex>(griddedData);
367 4 : lattice=arrayLattice;
368 : //logIO() << LogIO::DEBUGGING
369 : // << "Starting grid correction and FFT of image" << LogIO::POST;
370 :
371 : // Do the Grid-correction.
372 : {
373 4 : Vector<Complex> correction(nx);
374 4 : correction=Complex(1.0, 0.0);
375 : // Do the Grid-correction
376 4 : IPosition cursorShape(4, nx, 1, 1, 1);
377 4 : IPosition axisPath(4, 0, 1, 2, 3);
378 4 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
379 4 : LatticeIterator<Complex> lix(*lattice, lsx);
380 10244 : for(lix.reset();!lix.atEnd();lix++) {
381 10240 : gridder->correctX1D(correction, lix.position()(1));
382 10240 : lix.rwVectorCursor()/=correction;
383 : }
384 4 : }
385 4 : image->clearCache();
386 : // Now do the FFT2D in place
387 4 : LatticeFFT::cfft2d(*lattice);
388 :
389 : //logIO() << LogIO::DEBUGGING
390 : // << "Finished grid correction and FFT of image" << LogIO::POST;
391 :
392 4 : }
393 :
394 :
395 4 : void GridFT::finalizeToVis()
396 : {
397 :
398 : //cerr <<"Time to degrid " << timedegrid_p << endl;
399 4 : timedegrid_p=0.0;
400 :
401 4 : if(!arrayLattice.null()) arrayLattice=0;
402 4 : if(!lattice.null()) lattice=0;
403 4 : griddedData.resize();
404 :
405 4 : 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 4 : }
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 1 : 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 1 : }
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 0 : }
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 200 : Cube<Complex> data;
719 : //Fortran gridder need the flag as ints
720 200 : Cube<Int> flags;
721 200 : 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 200 : 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 200 : Matrix<Double> uvw(3, vb.uvw().nelements());
764 200 : uvw=0.0;
765 200 : Vector<Double> dphase(vb.uvw().nelements());
766 200 : Cube<Int> loc(2, nvischan, nRow);
767 200 : Cube<Int> off(2, nvischan, nRow);
768 200 : 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 200 : 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 200 : ixsub=2;
853 200 : iysub=2;
854 : }
855 0 : 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 200 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
867 1000 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
868 800 : sumwgt[icounter].resize(sumWeight.shape());
869 800 : 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 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 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 200 : }
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 72 : void GridFT::get(VisBuffer& vb, Int row)
1003 : {
1004 :
1005 72 : gridOk(convSupport_p);
1006 : // If row is -1 then we pass through all rows
1007 : Int startRow, endRow, nRow;
1008 72 : if (row < 0) {
1009 72 : nRow=vb.nRow();
1010 72 : startRow=0;
1011 72 : 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 72 : Matrix<Double> uvw(3, vb.uvw().nelements());
1024 72 : uvw=0.0;
1025 72 : Vector<Double> dphase(vb.uvw().nelements());
1026 72 : dphase=0.0;
1027 : //NEGATING to correct for an image inversion problem
1028 25344 : for (Int i=startRow;i<=endRow;i++) {
1029 75816 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
1030 25272 : uvw(2,i)=vb.uvw()(i)(2);
1031 : }
1032 72 : rotateUVW(uvw, dphase, vb);
1033 72 : 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 72 : 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 72 : if(doConversion_p[vb.spectralWindow()]){
1046 0 : matchChannel(vb.spectralWindow(), vb);
1047 : }
1048 : else{
1049 72 : chanMap.resize();
1050 72 : 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 72 : if(max(chanMap)==-1)
1056 0 : return;
1057 :
1058 72 : Cube<Complex> data;
1059 72 : Cube<Int> flags;
1060 72 : 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 72 : Int nvp=data.shape()(0);
1066 72 : Int nvc=data.shape()(1);
1067 72 : Int nvisrow=data.shape()(2);
1068 :
1069 :
1070 : //cerr << "get flags " << min(flags) << " " << max(flags) << endl;
1071 : Complex *datStorage;
1072 : Bool isCopy;
1073 72 : datStorage=data.getStorage(isCopy);
1074 :
1075 : ///
1076 72 : Cube<Int> loc(2, nvc, nvisrow);
1077 72 : Cube<Int> off(2, nvc, nRow);
1078 72 : Matrix<Complex> phasor(nvc, nRow);
1079 72 : Int csamp=convSampling_p;
1080 : Bool delphase;
1081 : Bool del;
1082 72 : Complex * phasorstor=phasor.getStorage(delphase);
1083 72 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1084 72 : const Double * scalestor=uvScale.getStorage(del);
1085 72 : const Double * offsetstor=uvOffset.getStorage(del);
1086 72 : const Double* uvstor= uvw.getStorage(del);
1087 72 : Int * locstor=loc.getStorage(del);
1088 72 : Int * offstor=off.getStorage(del);
1089 72 : const Double *dpstor=dphase.getStorage(del);
1090 : //Vector<Double> f1=interpVisFreq_p.copy();
1091 72 : Int nvchan=nvc;
1092 : Int irow;
1093 72 : Double cinv=Double(1.0)/C::c;
1094 72 : Int dow=0;
1095 72 : Int nth=1;
1096 : #ifdef _OPENMP
1097 72 : if(numthreads_p >0){
1098 0 : nth=min(numthreads_p, omp_get_max_threads());
1099 : }
1100 : else{
1101 72 : nth= omp_get_max_threads();
1102 : }
1103 72 : nth=min(4,nth);
1104 : #endif
1105 :
1106 :
1107 72 : #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 72 : Int rbeg=startRow+1;
1122 72 : Int rend=endRow+1;
1123 :
1124 :
1125 72 : Vector<Int> rowFlags(vb.nRow());
1126 72 : rowFlags=0;
1127 72 : rowFlags(vb.flagRow())=true;
1128 : //cerr << "rowFlags " << rowFlags << endl;
1129 72 : if(!usezero_p) {
1130 25344 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1131 25272 : 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 72 : const Complex* gridstor=griddedData.getStorage(delgrid);
1141 72 : const Double * convfuncstor=convFunc_p.getStorage(del);;
1142 :
1143 72 : const Int* pmapstor=polMap.getStorage(del);
1144 72 : const Int *cmapstor=chanMap.getStorage(del);
1145 72 : Int nc=nchan;
1146 72 : Int np=npol;
1147 72 : Int nxp=nx;
1148 72 : Int nyp=ny;
1149 72 : Int csupp=convSupport_p;
1150 72 : const Int * flagstor=flags.getStorage(del);
1151 72 : const Int * rowflagstor=rowFlags.getStorage(del);
1152 :
1153 :
1154 72 : Int npart=1;
1155 72 : if (nth >3){
1156 72 : npart=4;
1157 : }
1158 0 : else if(nth >1){
1159 0 : npart=2;
1160 : }
1161 :
1162 72 : Int ix=0;
1163 72 : #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 72 : data.putStorage(datStorage, isCopy);
1191 72 : griddedData.freeStorage(gridstor, delgrid);
1192 : //cerr << "Get min max " << min(data) << " " << max(data) << endl;
1193 : }
1194 72 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1195 :
1196 72 : }
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 0 : }
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 1 : Vector<Complex> correction(inx);
1271 1 : correction=Complex(1.0, 0.0);
1272 : // Do the Grid-correction
1273 1 : IPosition cursorShape(4, inx, 1, 1, 1);
1274 1 : IPosition axisPath(4, 0, 1, 2, 3);
1275 1 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1276 1 : 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 1 : }
1297 :
1298 1 : if(!isTiled) {
1299 : // Check the section from the image BEFORE converting to a lattice
1300 2 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1301 1 : IPosition stride(4, 1);
1302 2 : IPosition trc(blc+image->shape()-stride);
1303 : // Do the copy
1304 1 : IPosition start(4, 0);
1305 1 : image->put(griddedData(blc, trc));
1306 1 : }
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 0 : }
1382 :
1383 0 : Bool GridFT::fromRecord(String& error,
1384 : const RecordInterface& inRec)
1385 : {
1386 0 : Bool retval = true;
1387 0 : if(!FTMachine::fromRecord(error, inRec))
1388 0 : return false;
1389 0 : gridder=0; imageCache=0; lattice=0; arrayLattice=0;
1390 : //For some reason int64 seems to be getting lost...
1391 : //this is not an important parameter...so filing a jira
1392 0 : if(inRec.isDefined("cachesize"))
1393 0 : cachesize=inRec.asInt64("cachesize");
1394 : else
1395 0 : cachesize=1000000;
1396 0 : inRec.get("tilesize", tilesize);
1397 0 : inRec.get("convtype", convType);
1398 0 : inRec.get("uvscale", uvScale);
1399 0 : inRec.get("uvoffset", uvOffset);
1400 0 : inRec.get("istiled", isTiled);
1401 0 : inRec.get("maxabsdata", maxAbsData);
1402 0 : Vector<Int> center_loc(4), offset_loc(4);
1403 0 : inRec.get("centerloc", center_loc);
1404 0 : inRec.get("offsetloc", offset_loc);
1405 0 : uInt ndim4 = 4;
1406 0 : centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2),
1407 0 : center_loc(3));
1408 0 : offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2),
1409 0 : offset_loc(3));
1410 0 : inRec.get("padding", padding_p);
1411 0 : inRec.get("usezero", usezero_p);
1412 0 : inRec.get("nopadding", noPadding_p);
1413 :
1414 0 : machineName_p="GridFT";
1415 : ///setup some of the parameters
1416 0 : init();
1417 : /*if(!cmplxImage_p.null()){
1418 : //FTMachine::fromRecord would have recovered the image
1419 : // Might be changing the shape of sumWeight
1420 :
1421 : if(isTiled) {
1422 : lattice=CountedPtr<Lattice<Complex> >(image, false);
1423 : }
1424 : else {
1425 : // Make the grid the correct shape and turn it into an array lattice
1426 : // Check the section from the image BEFORE converting to a lattice
1427 : if(!toVis_p){
1428 : IPosition gridShape(4, nx, ny, npol, nchan);
1429 : griddedData.resize(gridShape);
1430 : griddedData=Complex(0.0);
1431 : }
1432 : else{
1433 : prepGridForDegrid();
1434 : }
1435 : }
1436 : };*/
1437 0 : return retval;
1438 0 : }
1439 :
1440 553 : void GridFT::ok() {
1441 553 : AlwaysAssert(image, AipsError);
1442 553 : }
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 1380 : String GridFT::name() const {
1514 :
1515 1380 : return machineName_p;
1516 :
1517 :
1518 : }
1519 :
1520 : } //# NAMESPACE CASA - END
1521 :
|
