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