Line data Source code
1 : //# WProjectFT.cc: Implementation of WProjectFT class
2 : //# Copyright (C) 2001-2016
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 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 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 :
29 : #include <casacore/casa/Quanta/UnitMap.h>
30 : #include <casacore/casa/Quanta/MVTime.h>
31 : #include <casacore/casa/Quanta/UnitVal.h>
32 : #include <casacore/casa/Utilities/CountedPtr.h>
33 : #include <casacore/measures/Measures/Stokes.h>
34 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
35 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
36 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
37 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
38 : #include <casacore/coordinates/Coordinates/Projection.h>
39 : #include <casacore/ms/MeasurementSets/MSColumns.h>
40 : #include <casacore/casa/BasicSL/Constants.h>
41 : #include <casacore/scimath/Mathematics/FFTServer.h>
42 : #include <synthesis/TransformMachines/WProjectFT.h>
43 : #include <synthesis/TransformMachines/WPConvFunc.h>
44 : #include <casacore/scimath/Mathematics/RigidVector.h>
45 : #include <msvis/MSVis/StokesVector.h>
46 : #include <synthesis/TransformMachines/StokesImageUtil.h>
47 : #include <msvis/MSVis/VisBuffer.h>
48 : #include <msvis/MSVis/VisSet.h>
49 : #include <msvis/MSVis/VisibilityIterator.h>
50 : #include <casacore/images/Images/ImageInterface.h>
51 : #include <casacore/images/Images/PagedImage.h>
52 : #include <casacore/casa/Containers/Block.h>
53 : #include <casacore/casa/Containers/Record.h>
54 : #include <casacore/casa/Arrays/ArrayLogical.h>
55 : #include <casacore/casa/Arrays/ArrayMath.h>
56 : #include <casacore/casa/Arrays/Array.h>
57 : #include <casacore/casa/Arrays/MaskedArray.h>
58 : #include <casacore/casa/Arrays/Vector.h>
59 : #include <casacore/casa/Arrays/Slice.h>
60 : #include <casacore/casa/Arrays/Matrix.h>
61 : #include <casacore/casa/Arrays/Cube.h>
62 : #include <casacore/casa/Arrays/MatrixIter.h>
63 : #include <casacore/casa/BasicSL/String.h>
64 : #include <casacore/casa/Utilities/Assert.h>
65 : #include <casacore/casa/Exceptions/Error.h>
66 : #include <iostream>
67 : #include <iomanip>
68 : #include <casacore/lattices/Lattices/ArrayLattice.h>
69 : #include <casacore/lattices/Lattices/SubLattice.h>
70 : #include <casacore/lattices/LRegions/LCBox.h>
71 : #include <casacore/lattices/LEL/LatticeExpr.h>
72 : #include <casacore/lattices/Lattices/LatticeCache.h>
73 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
74 : #include <casacore/lattices/Lattices/LatticeIterator.h>
75 : #include <casacore/lattices/Lattices/LatticeStepper.h>
76 : #include <casacore/casa/Utilities/CompositeNumber.h>
77 : #include <casacore/casa/OS/Timer.h>
78 : #include <casacore/casa/OS/HostInfo.h>
79 : #include <sstream>
80 : #ifdef _OPENMP
81 : #include <omp.h>
82 : #endif
83 : using namespace casacore;
84 : namespace casa { //# NAMESPACE CASA - BEGIN
85 :
86 0 : WProjectFT::WProjectFT( Int nWPlanes, Long icachesize, Int itilesize,
87 0 : Bool usezero, Bool useDoublePrec, const Double minW, const Double maxW, const Double rmsW)
88 : : FTMachine(), padding_p(1.0), nWPlanes_p(nWPlanes),
89 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
90 : gridder(0), isTiled(false),
91 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)), usezero_p(usezero),
92 0 : machineName_p("WProjectFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), minW_p(minW), maxW_p(maxW), rmsW_p(rmsW)
93 : {
94 0 : convSize=0;
95 0 : tangentSpecified_p=false;
96 0 : lastIndex_p=0;
97 0 : useDoubleGrid_p=useDoublePrec;
98 0 : wpConvFunc_p=new WPConvFunc(minW, maxW, rmsW);
99 :
100 0 : }
101 0 : WProjectFT::WProjectFT(Int nWPlanes,
102 : MPosition mLocation,
103 : Long icachesize, Int itilesize,
104 0 : Bool usezero, Float padding, Bool useDoublePrec, const Double minW, const Double maxW, const Double rmsW)
105 : : FTMachine(), padding_p(padding), nWPlanes_p(nWPlanes),
106 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
107 : gridder(0), isTiled(false),
108 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
109 : usezero_p(usezero),
110 0 : machineName_p("WProjectFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), minW_p(minW), maxW_p(maxW), rmsW_p(rmsW)
111 : {
112 0 : convSize=0;
113 0 : savedWScale_p=0.0;
114 0 : tangentSpecified_p=false;
115 0 : mLocation_p=mLocation;
116 0 : lastIndex_p=0;
117 0 : wpConvFunc_p=new WPConvFunc(minW, maxW, rmsW);
118 0 : useDoubleGrid_p=useDoublePrec;
119 0 : }
120 0 : WProjectFT::WProjectFT(
121 : Int nWPlanes, MDirection mTangent,
122 : MPosition mLocation,
123 : Long icachesize, Int itilesize,
124 0 : Bool usezero, Float padding, Bool useDoublePrec, const Double minW, const Double maxW, const Double rmsW)
125 : : FTMachine(), padding_p(padding), nWPlanes_p(nWPlanes),
126 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
127 : gridder(0), isTiled(false),
128 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
129 : usezero_p(usezero),
130 0 : machineName_p("WProjectFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0), minW_p(minW), maxW_p(maxW), rmsW_p(rmsW)
131 : {
132 0 : convSize=0;
133 0 : savedWScale_p=0.0;
134 0 : mTangent_p=mTangent;
135 0 : tangentSpecified_p=true;
136 0 : mLocation_p=mLocation;
137 0 : lastIndex_p=0;
138 0 : wpConvFunc_p=new WPConvFunc(minW, maxW, rmsW);
139 0 : useDoubleGrid_p=useDoublePrec;
140 0 : }
141 :
142 0 : WProjectFT::WProjectFT(const RecordInterface& stateRec)
143 0 : : FTMachine(),machineName_p("WProjectFT"), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0)
144 : {
145 : // Construct from the input state record
146 0 : String error;
147 0 : if (!fromRecord(error, stateRec)) {
148 0 : throw (AipsError("Failed to create WProjectFT: " + error));
149 : };
150 0 : }
151 : //----------------------------------------------------------------------
152 0 : WProjectFT& WProjectFT::operator=(const WProjectFT& other)
153 : {
154 0 : if(this!=&other) {
155 : //Do the base parameters
156 0 : FTMachine::operator=(other);
157 :
158 :
159 0 : padding_p=other.padding_p;
160 0 : nWPlanes_p=other.nWPlanes_p;
161 0 : imageCache=other.imageCache;
162 0 : cachesize=other.cachesize;
163 0 : tilesize=other.tilesize;
164 0 : if(other.gridder==0)
165 0 : gridder=0;
166 : else{
167 0 : uvScale.resize();
168 0 : uvOffset.resize();
169 0 : uvScale=other.uvScale;
170 0 : uvOffset=other.uvOffset;
171 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
172 0 : uvScale, uvOffset,
173 0 : "SF");
174 : }
175 :
176 0 : isTiled=other.isTiled;
177 : //lattice=other.lattice;
178 : //arrayLattice=other.arrayLattice;
179 0 : lattice=0;
180 0 : arrayLattice=0;
181 :
182 0 : maxAbsData=other.maxAbsData;
183 0 : centerLoc=other.centerLoc;
184 0 : offsetLoc=other.offsetLoc;
185 0 : usezero_p=other.usezero_p;
186 0 : machineName_p=other.machineName_p;
187 0 : wpConvFunc_p=other.wpConvFunc_p;
188 0 : timemass_p=0.0;
189 0 : timegrid_p=0.0;
190 0 : timedegrid_p=0.0;
191 0 : minW_p=other.minW_p;
192 0 : maxW_p=other.maxW_p;
193 0 : rmsW_p=other.rmsW_p;
194 : };
195 0 : return *this;
196 : };
197 :
198 : //----------------------------------------------------------------------
199 0 : WProjectFT::WProjectFT(const WProjectFT& other) :FTMachine(), machineName_p("WProjectFT"),
200 0 : timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0)
201 : {
202 0 : operator=(other);
203 0 : }
204 :
205 0 : FTMachine* WProjectFT::cloneFTM(){
206 0 : return new WProjectFT(*this);
207 : }
208 :
209 : //----------------------------------------------------------------------
210 0 : void WProjectFT::init() {
211 : /* if((padding_p*padding_p*image->shape().product())>cachesize) {
212 : isTiled=true;
213 : nx = image->shape()(0);
214 : ny = image->shape()(1);
215 : npol = image->shape()(2);
216 : nchan = image->shape()(3);
217 : }
218 : else {*/
219 : // We are padding.
220 0 : isTiled=false;
221 0 : CompositeNumber cn(uInt(image->shape()(0)*2));
222 0 : nx = cn.nextLargerEven(Int(padding_p*Float(image->shape()(0))-0.5));
223 0 : ny = cn.nextLargerEven(Int(padding_p*Float(image->shape()(1))-0.5));
224 : //cerr << "INIT shape " << image->shape()(0) << " " << image->shape()(1) << " new " << nx << " " << ny << endl;
225 0 : npol = image->shape()(2);
226 0 : nchan = image->shape()(3);
227 : //}
228 :
229 : // if(image->shape().product()>cachesize) {
230 : // isTiled=true;
231 : // }
232 : // else {
233 : // isTiled=false;
234 : // }
235 : //The Tiled version need some fixing: sof or now
236 0 : isTiled=false;
237 :
238 :
239 0 : sumWeight.resize(npol, nchan);
240 :
241 0 : wConvSize=max(0, nWPlanes_p);
242 0 : convSupport.resize(wConvSize);
243 0 : convSupport=0;
244 :
245 0 : uvScale.resize(3);
246 0 : uvScale=0.0;
247 0 : uvScale(0)=Float(nx)*image->coordinates().increment()(0);
248 0 : uvScale(1)=Float(ny)*image->coordinates().increment()(1);
249 0 : if(savedWScale_p==0.0){
250 0 : uvScale(2)=Float(wConvSize)*abs(image->coordinates().increment()(0));
251 : }
252 : else{
253 0 : uvScale(2)=savedWScale_p;
254 : }
255 0 : uvOffset.resize(3);
256 0 : uvOffset(0)=nx/2;
257 0 : uvOffset(1)=ny/2;
258 0 : uvOffset(2)=0;
259 :
260 :
261 :
262 0 : if(gridder) delete gridder; gridder=0;
263 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
264 0 : uvScale, uvOffset,
265 0 : "SF");
266 :
267 : // Set up image cache needed for gridding.
268 0 : if(imageCache) delete imageCache; imageCache=0;
269 :
270 : // The tile size should be large enough that the
271 : // extended convolution function can fit easily
272 0 : if(isTiled) {
273 0 : Float tileOverlap=0.5;
274 0 : tilesize=min(256,tilesize);
275 0 : IPosition tileShape=IPosition(4,tilesize,tilesize,npol,nchan);
276 0 : Vector<Float> tileOverlapVec(4);
277 0 : tileOverlapVec=0.0;
278 0 : tileOverlapVec(0)=tileOverlap;
279 0 : tileOverlapVec(1)=tileOverlap;
280 0 : Int tmpCacheVal=static_cast<Int>(cachesize);
281 0 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
282 : tileOverlapVec,
283 0 : (tileOverlap>0.0));
284 :
285 : }
286 0 : }
287 :
288 : // This is nasty, we should use CountedPointers here.
289 0 : WProjectFT::~WProjectFT() {
290 0 : if(imageCache) delete imageCache; imageCache=0;
291 0 : if(gridder) delete gridder; gridder=0;
292 : /*
293 : if(arrayLattice) delete arrayLattice; arrayLattice=0;
294 :
295 : Int numofmodels=convFunctions_p.nelements();
296 : for (Int k=0; k< numofmodels; ++k){
297 : delete convFunctions_p[k];
298 : delete convSupportBlock_p[k];
299 :
300 : }
301 : */
302 : // convFuctions_p.resize();
303 : // convSupportBlock_p.resize();
304 :
305 0 : }
306 :
307 :
308 0 : void WProjectFT::setConvFunc(CountedPtr<WPConvFunc>& pbconvFunc){
309 :
310 0 : wpConvFunc_p=pbconvFunc;
311 0 : }
312 0 : CountedPtr<WPConvFunc>& WProjectFT::getConvFunc(){
313 :
314 0 : return wpConvFunc_p;
315 : }
316 :
317 0 : void WProjectFT::findConvFunction(const ImageInterface<Complex>& image,
318 : const VisBuffer& vb) {
319 :
320 :
321 :
322 :
323 :
324 0 : wpConvFunc_p->findConvFunction(image, vb, wConvSize, uvScale, uvOffset,
325 0 : padding_p,
326 0 : convSampling,
327 0 : convFunc, convSize, convSupport,
328 0 : savedWScale_p);
329 :
330 0 : nWPlanes_p=convSupport.nelements();
331 0 : wConvSize=max(1,nWPlanes_p);
332 0 : uvScale(2)=savedWScale_p;
333 :
334 :
335 :
336 0 : }
337 :
338 0 : void WProjectFT::initializeToVis(ImageInterface<Complex>& iimage,
339 : const VisBuffer& vb)
340 : {
341 0 : image=&iimage;
342 0 : toVis_p=true;
343 0 : ok();
344 :
345 : // if(convSize==0) {
346 0 : init();
347 : // }
348 0 : findConvFunction(*image, vb);
349 :
350 :
351 : // Initialize the maps for polarization and channel. These maps
352 : // translate visibility indices into image indices
353 0 : initMaps(vb);
354 :
355 : // nx = image->shape()(0);
356 : // ny = image->shape()(1);
357 : // npol = image->shape()(2);
358 : // nchan = image->shape()(3);
359 :
360 :
361 0 : isTiled=false;
362 : // If we are memory-based then read the image in and create an
363 : // ArrayLattice otherwise just use the PagedImage
364 : /*if(isTiled) {
365 : lattice=CountedPtr<Lattice<Complex> > (image, false);
366 : }
367 : else {
368 : }
369 : */
370 : //AlwaysAssert(lattice, AipsError);
371 :
372 0 : prepGridForDegrid();
373 0 : }
374 :
375 0 : void WProjectFT::prepGridForDegrid(){
376 0 : IPosition gridShape(4, nx, ny, npol, nchan);
377 0 : griddedData.resize(gridShape);
378 0 : griddedData=Complex(0.0);
379 :
380 :
381 0 : IPosition stride(4, 1);
382 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
383 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
384 0 : IPosition trc(blc+image->shape()-stride);
385 :
386 0 : IPosition start(4, 0);
387 0 : griddedData(blc, trc) = image->getSlice(start, image->shape());
388 :
389 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
390 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
391 0 : lattice=arrayLattice;
392 :
393 0 : logIO() << LogIO::DEBUGGING << "Starting FFT of image" << LogIO::POST;
394 :
395 :
396 0 : Vector<Float> sincConvX(nx);
397 0 : for (Int ix=0;ix<nx;ix++) {
398 0 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
399 0 : if(ix==nx/2) {
400 0 : sincConvX(ix)=1.0;
401 : }
402 : else {
403 0 : sincConvX(ix)=sin(x)/x;
404 : }
405 : }
406 0 : Vector<Float> sincConvY(ny);
407 0 : for (Int ix=0;ix<ny;ix++) {
408 0 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
409 0 : if(ix==ny/2) {
410 0 : sincConvY(ix)=1.0;
411 : }
412 : else {
413 0 : sincConvY(ix)=sin(x)/x;
414 : }
415 : }
416 :
417 :
418 0 : Vector<Complex> correction(nx);
419 0 : correction=Complex(1.0, 0.0);
420 : // Do the Grid-correction
421 0 : IPosition cursorShape(4, nx, 1, 1, 1);
422 0 : IPosition axisPath(4, 0, 1, 2, 3);
423 0 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
424 0 : LatticeIterator<Complex> lix(*lattice, lsx);
425 0 : for(lix.reset();!lix.atEnd();lix++) {
426 0 : Int iy=lix.position()(1);
427 0 : gridder->correctX1D(correction,iy);
428 0 : for (Int ix=0;ix<nx;ix++) {
429 0 : correction(ix)/=(sincConvX(ix)*sincConvY(iy));
430 : }
431 0 : lix.rwVectorCursor()/=correction;
432 : }
433 :
434 : // Now do the FFT2D in place
435 0 : LatticeFFT::cfft2d(*lattice);
436 :
437 0 : logIO() << LogIO::DEBUGGING << "Finished FFT" << LogIO::POST;
438 :
439 :
440 :
441 0 : }
442 :
443 :
444 0 : void WProjectFT::finalizeToVis()
445 : {
446 :
447 : //cerr <<"Time to degrid " << timedegrid_p << endl;
448 0 : timedegrid_p=0.0;
449 0 : if(!arrayLattice.null()) arrayLattice=0;
450 0 : if(!lattice.null()) lattice=0;
451 0 : griddedData.resize();
452 0 : if(isTiled) {
453 :
454 0 : logIO() << LogOrigin("WProjectFT", "finalizeToVis") << LogIO::NORMAL;
455 :
456 0 : AlwaysAssert(imageCache, AipsError);
457 0 : AlwaysAssert(image, AipsError);
458 0 : ostringstream o;
459 0 : imageCache->flush();
460 0 : imageCache->showCacheStatistics(o);
461 0 : logIO() << o.str() << LogIO::POST;
462 : }
463 0 : }
464 :
465 :
466 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
467 : // grid.
468 0 : void WProjectFT::initializeToSky(ImageInterface<Complex>& iimage,
469 : Matrix<Float>& weight,
470 : const VisBuffer& vb)
471 : {
472 : // image always points to the image
473 0 : image=&iimage;
474 0 : toVis_p=false;
475 :
476 : // if(convSize==0) {
477 0 : init();
478 : // }
479 0 : findConvFunction(*image, vb);
480 :
481 :
482 : // Initialize the maps for polarization and channel. These maps
483 : // translate visibility indices into image indices
484 0 : initMaps(vb);
485 :
486 : // nx = image->shape()(0);
487 : // ny = image->shape()(1);
488 : // npol = image->shape()(2);
489 : // nchan = image->shape()(3);
490 :
491 : // if(image->shape().product()>cachesize) {
492 : // isTiled=true;
493 : // }
494 : // else {
495 : // isTiled=false;
496 : // }
497 0 : isTiled=false;
498 0 : sumWeight=0.0;
499 0 : weight.resize(sumWeight.shape());
500 0 : weight=0.0;
501 :
502 : // Initialize for in memory or to disk gridding. lattice will
503 : // point to the appropriate Lattice, either the ArrayLattice for
504 : // in memory gridding or to the image for to disk gridding.
505 0 : if(isTiled) {
506 0 : imageCache->flush();
507 0 : image->set(Complex(0.0));
508 0 : lattice=CountedPtr<Lattice<Complex> > (image, false);
509 : }
510 : else {
511 0 : IPosition gridShape(4, nx, ny, npol, nchan);
512 0 : if( !useDoubleGrid_p )
513 : {
514 0 : griddedData.resize(gridShape);
515 0 : griddedData=Complex(0.0);
516 0 : griddedData2.resize();
517 : }
518 : else
519 0 : {griddedData.resize();
520 0 : griddedData2.resize(gridShape);
521 0 : griddedData2=DComplex(0.0);
522 : }
523 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
524 :
525 : }
526 : //AlwaysAssert(lattice, AipsError);
527 :
528 0 : }
529 :
530 0 : void WProjectFT::finalizeToSky()
531 : {
532 :
533 :
534 : //cerr <<"Time to massage data " << timemass_p << endl;
535 : //cerr <<"Time to grid data " << timegrid_p << endl;
536 0 : timemass_p=0.0;
537 0 : timegrid_p=0.0;
538 : // Now we flush the cache and report statistics
539 : // For memory based, we don't write anything out yet.
540 0 : if(isTiled) {
541 0 : logIO() << LogOrigin("WProjectFT", "finalizeToSky") << LogIO::NORMAL;
542 :
543 0 : AlwaysAssert(image, AipsError);
544 0 : AlwaysAssert(imageCache, AipsError);
545 0 : imageCache->flush();
546 0 : ostringstream o;
547 0 : imageCache->showCacheStatistics(o);
548 0 : logIO() << o.str() << LogIO::POST;
549 : }
550 0 : }
551 :
552 0 : Array<Complex>* WProjectFT::getDataPointer(const IPosition& centerLoc2D,
553 : Bool readonly) {
554 : Array<Complex>* result;
555 : // Is tiled: get tiles and set up offsets
556 0 : centerLoc(0)=centerLoc2D(0);
557 0 : centerLoc(1)=centerLoc2D(1);
558 0 : result=&imageCache->tile(offsetLoc, centerLoc, readonly);
559 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
560 0 : return result;
561 : }
562 :
563 : #define NEED_UNDERSCORES
564 : #if defined(NEED_UNDERSCORES)
565 : #define gwgrid gwgrid_
566 : #define gwproj gwproj_
567 : #define dwproj dwproj_
568 : #define sectgwgridd sectgwgridd_
569 : #define sectgwgrids sectgwgrids_
570 : #define sectdwgrid sectdwgrid_
571 : #define locuvw locuvw_
572 : #endif
573 :
574 : extern "C" {
575 : void locuvw(const Double*, const Double*, const Double*, const Int*, const Double*, const Double*, const Int*,
576 : Int*, Int*, Complex*, const Int*, const Int*, const Double*);
577 : //Double precision gridding
578 : void gwgrid(const Double*,
579 : Double*,
580 : const Complex*,
581 : Int*,
582 : Int*,
583 : Int*,
584 : const Int*,
585 : const Int*,
586 : const Float*,
587 : Int*,
588 : Int*,
589 : Double*,
590 : Double*,
591 : DComplex*,
592 : Int*,
593 : Int*,
594 : Int *,
595 : Int *,
596 : const Double*,
597 : const Double*,
598 : Int*,
599 : Int*,
600 : Int*,
601 : Int*,
602 : const Complex*,
603 : Int*,
604 : Int*,
605 : Double*);
606 :
607 : void sectgwgridd(const Double*,
608 : const Complex*,
609 : const Int*,
610 : const Int*,
611 : const Int*,
612 : const Int*,
613 : const Int*,
614 : const Float*,
615 : const Int*,
616 : DComplex*,
617 : const Int*,
618 : const Int*,
619 : const Int *,
620 : const Int *,
621 : //support
622 : const Int*,
623 : const Int*,
624 : const Int*,
625 : const Int*,
626 : const Complex*,
627 : const Int*,
628 : const Int*,
629 : Double*,
630 : //x0
631 : const Int*,
632 : const Int*,
633 : const Int*,
634 : const Int*,
635 : const Int*,
636 : const Int*,
637 : const Int*,
638 : const Int*,
639 : const Complex*
640 : );
641 :
642 : //Single precision gridding
643 : void sectgwgrids(const Double*,
644 : const Complex*,
645 : const Int*,
646 : const Int*,
647 : const Int*,
648 : const Int*,
649 : const Int*,
650 : const Float*,
651 : const Int*,
652 : Complex*,
653 : const Int*,
654 : const Int*,
655 : const Int *,
656 : const Int *,
657 : //support
658 : const Int*,
659 : const Int*,
660 : const Int*,
661 : const Int*,
662 : const Complex*,
663 : const Int*,
664 : const Int*,
665 : Double*,
666 : //x0
667 : const Int*,
668 : const Int*,
669 : const Int*,
670 : const Int*,
671 : const Int*,
672 : const Int*,
673 : const Int*,
674 : const Int*,
675 : const Complex*
676 : );
677 :
678 :
679 : void gwproj(const Double*,
680 : Double*,
681 : const Complex*,
682 : Int*,
683 : Int*,
684 : Int*,
685 : const Int*,
686 : const Int*,
687 : const Float*,
688 : Int*,
689 : Int*,
690 : Double*,
691 : Double*,
692 : Complex*,
693 : Int*,
694 : Int*,
695 : Int *,
696 : Int *,
697 : const Double*,
698 : const Double*,
699 : Int*,
700 : Int*,
701 : Int*,
702 : Int*,
703 : const Complex*,
704 : Int*,
705 : Int*,
706 : Double*);
707 :
708 : void sectdwgrid(const Double*,
709 : Complex*,
710 : const Int*,
711 : const Int*,
712 : const Int*,
713 : const Int*,
714 : const Int*,
715 : const Complex*,
716 : const Int*,
717 : const Int*,
718 : const Int *,
719 : const Int *,
720 : //support
721 : const Int*,
722 : const Int*,
723 : const Int*,
724 : const Int*,
725 : const Complex*,
726 : const Int*,
727 : const Int*,
728 : //rbeg, rend, loc, off, phasor
729 : const Int*,
730 : const Int*,
731 : const Int*,
732 : const Int*,
733 : const Complex*);
734 : void dwproj(const Double*,
735 : Double*,
736 : Complex*,
737 : Int*,
738 : Int*,
739 : const Int*,
740 : const Int*,
741 : Int*,
742 : Int*,
743 : Double*,
744 : Double*,
745 : const Complex*,
746 : Int*,
747 : Int*,
748 : Int *,
749 : Int *,
750 : const Double*,
751 : const Double*,
752 : Int*,
753 : Int*,
754 : Int*,
755 : Int*,
756 : const Complex*,
757 : Int*,
758 : Int*);
759 : }
760 0 : void WProjectFT::put(const VisBuffer& vb, Int row, Bool dopsf,
761 : FTMachine::Type type)
762 : {
763 :
764 :
765 : //Check if ms has changed then cache new spw and chan selection
766 0 : if(vb.newMS())
767 0 : matchAllSpwChans(vb);
768 :
769 : //Here we redo the match or use previous match
770 :
771 : //Channel matching for the actual spectral window of buffer
772 0 : if(doConversion_p[vb.spectralWindow()]){
773 0 : matchChannel(vb.spectralWindow(), vb);
774 : }
775 : else{
776 0 : chanMap.resize();
777 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
778 : }
779 :
780 : //No point in reading data if its not matching in frequency
781 0 : if(max(chanMap)==-1)
782 0 : return;
783 0 : Timer tim;
784 0 : tim.mark();
785 :
786 : const Matrix<Float> *imagingweight;
787 0 : imagingweight=&(vb.imagingWeight());
788 :
789 0 : if(dopsf) type=FTMachine::PSF;
790 :
791 0 : Cube<Complex> data;
792 : //Fortran gridder need the flag as ints
793 0 : Cube<Int> flags;
794 0 : Matrix<Float> elWeight;
795 0 : interpolateFrequencyTogrid(vb, *imagingweight,data, flags, elWeight, type);
796 :
797 :
798 : Bool iswgtCopy;
799 : const Float *wgtStorage;
800 0 : wgtStorage=elWeight.getStorage(iswgtCopy);
801 :
802 :
803 : Bool isCopy;
804 0 : const Complex *datStorage=0;
805 :
806 0 : if(!dopsf)
807 0 : datStorage=data.getStorage(isCopy);
808 :
809 :
810 : // If row is -1 then we pass through all rows
811 : Int startRow, endRow, nRow;
812 0 : if (row==-1) {
813 0 : nRow=vb.nRow();
814 0 : startRow=0;
815 0 : endRow=nRow-1;
816 : } else {
817 0 : nRow=1;
818 0 : startRow=row;
819 0 : endRow=row;
820 : }
821 :
822 : // Get the uvws in a form that Fortran can use and do that
823 : // necessary phase rotation. On a Pentium Pro 200 MHz
824 : // when null, this step takes about 50us per uvw point. This
825 : // is just barely noticeable for Stokes I continuum and
826 : // irrelevant for other cases.
827 0 : Matrix<Double> uvw(3, vb.uvw().nelements());
828 0 : uvw=0.0;
829 0 : Vector<Double> dphase(vb.uvw().nelements());
830 0 : dphase=0.0;
831 : //NEGATING to correct for an image inversion problem
832 0 : for (Int i=startRow;i<=endRow;i++) {
833 0 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
834 0 : uvw(2,i)=vb.uvw()(i)(2);
835 : }
836 :
837 0 : rotateUVW(uvw, dphase, vb);
838 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
839 :
840 :
841 :
842 : // Take care of translation of Bools to Integer
843 0 : Int idopsf=0;
844 0 : if(dopsf) idopsf=1;
845 :
846 0 : Vector<Int> rowFlags(vb.nRow());
847 0 : rowFlags=0;
848 0 : rowFlags(vb.flagRow())=true;
849 0 : if(!usezero_p) {
850 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
851 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
852 : }
853 : }
854 :
855 :
856 : Bool del;
857 : // IPosition s(flags.shape());
858 0 : Vector<Int> s(flags.shape().nelements());
859 0 : convertArray(s, flags.shape().asVector());
860 0 : Int nvp=s[0];
861 0 : Int nvc=s[1];
862 0 : Int nvisrow=s[2];
863 :
864 0 : Int csamp=convSampling;
865 :
866 : Bool uvwcopy;
867 0 : const Double *uvwstor=uvw.getStorage(uvwcopy);
868 : Bool gridcopy;
869 : Bool convcopy;
870 0 : const Complex *convstor=convFunc.getStorage(convcopy);
871 0 : Cube<Int> loc(3, nvc, nRow);
872 0 : Cube<Int> off(3, nvc, nRow);
873 0 : Matrix<Complex> phasor(nvc, nRow);
874 : Bool delphase;
875 0 : Complex * phasorstor=phasor.getStorage(delphase);
876 0 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
877 0 : const Double * scalestor=uvScale.getStorage(del);
878 0 : const Double * offsetstor=uvOffset.getStorage(del);
879 : Bool delloc, deloff;
880 0 : Int * locstor=loc.getStorage(delloc);
881 0 : Int * offstor=off.getStorage(deloff);
882 0 : const Double *dpstor=dphase.getStorage(del);
883 0 : Double cinv=Double(1.0)/C::c;
884 : Int irow;
885 0 : Int dow=1;
886 0 : Int nth=1;
887 : #ifdef _OPENMP
888 0 : if(numthreads_p >0){
889 0 : nth=min(numthreads_p, omp_get_max_threads());
890 : }
891 : else{
892 0 : nth= omp_get_max_threads();
893 : }
894 : //nth=min(4,nth);
895 : #endif
896 :
897 :
898 :
899 :
900 0 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvc, scalestor, offsetstor, csamp, phasorstor, uvwstor, locstor, offstor, dpstor, cinv, dow) shared(startRow, endRow) num_threads(nth)
901 : {
902 : #pragma omp for schedule(dynamic)
903 : for (irow=startRow; irow<=endRow;irow++){
904 : //locateuvw(uvwstor,dpstor, visfreqstor, nvc, scalestor, offsetstor, csamp,
905 : // locstor,
906 : // offstor, phasorstor, irow, true);
907 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
908 : }
909 :
910 :
911 :
912 :
913 : }//end pragma parallel
914 : //Int maxx=0;
915 0 : Int minx=1000000;
916 : //Int maxy=0;
917 0 : Int miny=1000000;
918 : //Int maxsupp=max(convSupport);
919 : /*loc.putStorage(locstor, delloc);
920 : maxx=max(loc.yzPlane(0))+convSize;
921 : maxx=min(nx-1, maxx-1);
922 : minx=min(loc.yzPlane(0))-convSize;
923 : minx=max(0, minx-1);
924 : maxy=max(loc.yzPlane(1))+convSize;
925 : maxy=min(ny-1, maxy-1);
926 : miny=min(loc.yzPlane(1))-convSize;
927 : miny=max(0,miny-1);
928 : locstor=loc.getStorage(delloc);
929 : */
930 : //cerr << "dels " << delloc << " " << deloff << endl;
931 : //cerr << "LOC " << min(loc.xzPlane(0)) << " max " << loc.shape() << endl;
932 0 : timemass_p +=tim.real();
933 : Int ixsub, iysub, icounter;
934 0 : ixsub=1;
935 0 : iysub=1;
936 0 : if(nth > 4){
937 0 : ixsub=8;
938 0 : iysub=8;
939 : }
940 : else{
941 0 : ixsub=2;
942 0 : iysub=2;
943 : }
944 : /* else if(nth >1){
945 : ixsub=2;
946 : iysub=1;
947 : }
948 : */
949 :
950 : //nxsub=nx;
951 : //nysub=ny;
952 : //nxsub=maxx-minx+1;
953 : //nysub=maxy-miny+1;
954 0 : Int rbeg=startRow+1;
955 0 : Int rend=endRow+1;
956 0 : const Int* pmapstor=polMap.getStorage(del);
957 0 : const Int *cmapstor=chanMap.getStorage(del);
958 0 : Int nc=nchan;
959 0 : Int np=npol;
960 : ////For now no limitation
961 0 : minx=0;
962 0 : miny=0;
963 0 : Int nxp=nx;
964 0 : Int nyp=ny;
965 :
966 :
967 0 : Int nxcopy=nx;
968 0 : Int nycopy=ny;
969 :
970 0 : Int csize=convSize;
971 0 : Int wcsize=wConvSize;
972 0 : const Int * flagstor=flags.getStorage(del);
973 0 : const Int * rowflagstor=rowFlags.getStorage(del);
974 0 : const Int * suppstor=convSupport.getStorage(del);
975 : //for (Int biba=4; biba < 60; biba+=2){
976 : // cerr << "biba 0 " << biba << endl;
977 : // ixsub=biba;
978 : // iysub=biba;
979 :
980 0 : tim.mark();
981 : //clock_t tick;
982 : //tick = clock();
983 : //time_t now;
984 : //time(&now);
985 :
986 0 : Vector<Int> x0p, y0p, nxpsub, nypsub;
987 :
988 0 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
989 : // if(timegrid_p==0.0 && dopsf){
990 : {
991 0 : xsect_p.resize(ixsub*iysub);
992 0 : ysect_p.resize(ixsub*iysub);
993 0 : nxsect_p.resize(ixsub*iysub);
994 0 : nysect_p.resize(ixsub*iysub);
995 0 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
996 : ///testoo
997 : //minx=0;
998 : //miny=0;
999 : //nxp=nx;
1000 : //nyp=ny;
1001 : ////////
1002 0 : findGridSector(nxp, nyp, ixsub, iysub, minx, miny, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
1003 : }
1004 : }
1005 : /*if((Int(timegrid_p)%10000)==0){
1006 : cerr << timegrid_p << endl;
1007 : cerr << "xsect " << xsect_p << "\n ysect " << ysect_p << "\n nxsect "<< nxsect_p << "\n nysect " << nysect_p << endl;
1008 :
1009 : }
1010 : */
1011 0 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1012 0 : sumwgt[icounter].resize(sumWeight.shape());
1013 0 : sumwgt[icounter].set(0.0);
1014 : }
1015 0 : Vector<Int> xsect, ysect, nxsect, nysect;
1016 0 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
1017 : //cerr << "useDoubleGrid " << useDoubleGrid_p << endl;
1018 0 : if(!useDoubleGrid_p){
1019 0 : Complex *gridstor=griddedData.getStorage(gridcopy);
1020 0 : #pragma omp parallel default(none) private(icounter, del) firstprivate(idopsf, uvwstor, datStorage, wgtStorage, flagstor, rowflagstor, convstor, pmapstor, cmapstor, gridstor, suppstor, nxp, nyp, nxcopy, nycopy, np, nc,ixsub, iysub, rend, rbeg, csamp, csize, wcsize, nvp, nvc, nvisrow, phasorstor, locstor, offstor, minx, miny, xsect, ysect, nxsect, nysect) shared(sumwgt) num_threads(nth)
1021 : {
1022 : #pragma omp for schedule(dynamic)
1023 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
1024 : //Int realicounter=icounter%2==0 ? ixsub*iysub/2+icounter/2 : ixsub*iysub/2-icounter/2-1;
1025 : //findGridSector(nxp, nyp, ixsub, iysub, minx, miny, icounter, x0, y0, nxsub, nysub, true);
1026 : Int x0=xsect(icounter);
1027 : Int y0=ysect(icounter);
1028 : Int nxsub=nxsect(icounter);
1029 : Int nysub=nysect(icounter);
1030 :
1031 : sectgwgrids(uvwstor,
1032 : datStorage,
1033 : &nvp,
1034 : &nvc,
1035 : &idopsf,
1036 : flagstor,
1037 : rowflagstor,
1038 : wgtStorage,
1039 : &nvisrow,
1040 : gridstor,
1041 : &nxcopy,
1042 : &nycopy,
1043 : &np,
1044 : &nc,
1045 : suppstor,
1046 : &csize,
1047 : &csamp,
1048 : &wcsize,
1049 : convstor,
1050 : cmapstor,
1051 : pmapstor,
1052 : (sumwgt[icounter]).getStorage(del),
1053 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
1054 : phasorstor);
1055 : }
1056 : }//end pragma parallel
1057 :
1058 0 : timegrid_p+=tim.real();
1059 : //if(dopsf)
1060 : // tweakGridSector(nx, ny, ixsub, iysub, xsect_p, ysect_p, nxsect_p, nysect_p);
1061 0 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1062 : //cerr << "ixsub, iysub " << icounter%ixsub << " " << icounter/ixsub << " sumwgt " << sumwgt[icounter](0,0) << endl;
1063 0 : sumWeight=sumWeight+sumwgt[icounter];
1064 : }
1065 0 : griddedData.putStorage(gridstor, gridcopy);
1066 :
1067 : }
1068 : else{
1069 0 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
1070 0 : #pragma omp parallel default(none) private(icounter, del) firstprivate(idopsf, uvwstor, datStorage, wgtStorage, flagstor, rowflagstor, convstor, pmapstor, cmapstor, gridstor, suppstor, nxp, nyp, nxcopy, nycopy, np, nc,ixsub, iysub, rend, rbeg, csamp, csize, wcsize, nvp, nvc, nvisrow, phasorstor, locstor, offstor,minx,miny, xsect, ysect, nxsect, nysect) shared(sumwgt) num_threads(nth)
1071 : {
1072 : #pragma omp for schedule(dynamic)
1073 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
1074 : //Int realicounter=icounter%2==0 ? ixsub*iysub/2+icounter/2 : ixsub*iysub/2-icounter/2-1;
1075 : // findGridSector(nxp, nyp, ixsub, iysub, minx, miny, icounter, x0, y0, nxsub, nysub, true);
1076 : //Double wt0=omp_get_wtime();
1077 : Int x0=xsect(icounter);
1078 : Int y0=ysect(icounter);
1079 : Int nxsub=nxsect(icounter);
1080 : Int nysub=nysect(icounter);
1081 :
1082 : sectgwgridd(uvwstor,
1083 : datStorage,
1084 : &nvp,
1085 : &nvc,
1086 : &idopsf,
1087 : flagstor,
1088 : rowflagstor,
1089 : wgtStorage,
1090 : &nvisrow,
1091 : gridstor,
1092 : &nxcopy,
1093 : &nycopy,
1094 : &np,
1095 : &nc,
1096 : suppstor,
1097 : &csize,
1098 : &csamp,
1099 : &wcsize,
1100 : convstor,
1101 : cmapstor,
1102 : pmapstor,
1103 : (sumwgt[icounter]).getStorage(del),
1104 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
1105 : phasorstor);
1106 : //cerr << "icounter " << icounter << " " << icounter%ixsub << " " << icounter/ixsub << " " << omp_get_wtime()-wt0 << endl;
1107 : }
1108 : }//end pragma parallel
1109 0 : timegrid_p+=1.0;
1110 : //timegrid_p+=tim.real();
1111 : //if(dopsf)
1112 : // tweakGridSector(nx, ny, ixsub, iysub, xsect_p, ysect_p, nxsect_p, nysect_p);
1113 0 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1114 : //cerr << "ixsub, iysub " << icounter%ixsub << " " << icounter/ixsub << " sumwgt " << sumwgt[icounter](0,0) << endl;
1115 0 : sumWeight=sumWeight+sumwgt[icounter];
1116 : }
1117 0 : griddedData2.putStorage(gridstor, gridcopy);
1118 : }
1119 :
1120 :
1121 :
1122 : // std::cerr << std::setprecision(10) << "biba " << biba << " time " << tim.real() << " timall " << tim.all_usec() << endl;
1123 : // }
1124 0 : uvw.freeStorage(uvwstor, uvwcopy);
1125 0 : convFunc.freeStorage(convstor, convcopy);
1126 :
1127 0 : if(!dopsf)
1128 0 : data.freeStorage(datStorage, isCopy);
1129 0 : elWeight.freeStorage(wgtStorage,iswgtCopy);
1130 : //timegrid_p+=tim.real();
1131 : //timegrid_p+= Double(clock()-tick);
1132 : // time_t now1;
1133 : //time(&now1);
1134 : //timegrid_p+=difftime(now1, now);
1135 : }
1136 :
1137 0 : void WProjectFT::findGridSector(const Int& nxp, const Int& nyp, const Int& ixsub, const Int& iysub, const Int& minx, const Int& miny, const Int& icounter, Int& x0, Int& y0, Int& nxsub, Int& nysub, const Bool linear){
1138 : /* Vector<Int> ord(36);
1139 : ord(0)=14;
1140 : ord(1)=15;
1141 : ord(2)=20;
1142 : ord(3)=21;ord(4)=13;
1143 : ord(5)=16;ord(6)=19;ord(7)=22;ord(8)=8;ord(9)=9;
1144 : ord(10)=26;ord(11)=27;ord(12)=25;ord(13)=28;ord(14)=7;
1145 : ord(15)=10;ord(16)=32;ord(17)=33;ord(18)=2;ord(19)=3;
1146 : ord(20)=18;ord(21)=23;ord(22)=12;ord(23)=17;ord(24)=1;
1147 : ord(25)=4;ord(26)=6;ord(27)=11;ord(28)=24;ord(29)=29;
1148 : ord(30)=31;ord(31)=34;ord(32)=0;ord(33)=5;ord(34)=30;
1149 : ord(35)=35;
1150 : */
1151 : /*
1152 : Int ix= (icounter+1)-((icounter)/ixsub)*ixsub;
1153 : Int iy=(icounter)/ixsub+1;
1154 : y0=(nyp/iysub)*(iy-1)+1+miny;
1155 : nysub=nyp/iysub;
1156 : if( iy == iysub) {
1157 : nysub=nyp-(nyp/iysub)*(iy-1);
1158 : }
1159 : x0=(nxp/ixsub)*(ix-1)+1+minx;
1160 : nxsub=nxp/ixsub;
1161 : if( ix == ixsub){
1162 : nxsub=nxp-(nxp/ixsub)*(ix-1);
1163 : }
1164 : */
1165 :
1166 :
1167 : {
1168 0 : Int elrow=icounter/ixsub;
1169 0 : Int elcol=(icounter-elrow*ixsub);
1170 : //cerr << "row "<< elrow << " col " << elcol;
1171 : //nxsub=Int(floor(((ceil(fabs(float(2*elcol+1-ixsub)/2.0))-1.0)*5 +1)*nxp/36.0 + 0.5));
1172 0 : Float factor=0;
1173 0 : for (Int k=0; k < ixsub/2; ++k)
1174 0 : factor= linear ? factor+(k+1): factor+sqrt(Float(k+1));
1175 : //factor= linear ? factor+(k+1): factor+(k+1)*(k+1)*(k+1);
1176 0 : factor *= 2.0;
1177 0 : if(linear)
1178 0 : nxsub=Int(floor((ceil(fabs(float(2*elcol+1-ixsub)/2.0))/factor)*nxp + 0.5));
1179 : else
1180 : //nxsub=Int(floor((ceil(fabs(float(2*elcol+1-ixsub)/2.0))*ceil(fabs(float(2*elcol+1-ixsub)/2.0))*ceil(fabs(float(2*elcol+1-ixsub)/2.0))/factor)*nxp + 0.5));
1181 0 : nxsub=Int(floor((sqrt(ceil(fabs(float(2*elcol+1-ixsub)/2.0)))/factor)*nxp + 0.5));
1182 : //cerr << nxp << " col " << elcol << " nxsub " << nxsub << endl;
1183 0 : x0=minx;
1184 0 : elcol-=1;
1185 0 : while(elcol >= 0){
1186 : //x0+=Int(floor(((ceil(fabs(float(2*elcol+1-ixsub)/2.0))-1.0)*5 +1)*nxp/36.0+0.5));
1187 :
1188 0 : if(linear)
1189 0 : x0+=Int(floor((ceil(fabs(float(2*elcol+1-ixsub)/2.0))/factor)*nxp + 0.5));
1190 : else
1191 : //x0+=Int(floor((ceil(fabs(float(2*elcol+1-ixsub)/2.0))*ceil(fabs(float(2*elcol+1-ixsub)/2.0))*ceil(fabs(float(2*elcol+1-ixsub)/2.0))/factor)*nxp + 0.5));
1192 0 : x0+=Int(floor((sqrt(ceil(fabs(float(2*elcol+1-ixsub)/2.0)))/factor)*nxp + 0.5));
1193 0 : elcol-=1;
1194 : }
1195 0 : factor=0;
1196 0 : for (Int k=0; k < iysub/2; ++k)
1197 : //factor=linear ? factor+(k+1): factor+(k+1)*(k+1)*(k+1);
1198 0 : factor= linear ? factor+(k+1): factor+sqrt(Float(k+1));
1199 0 : factor *= 2.0;
1200 : //nysub=Int(floor(((ceil(fabs(float(2*elrow+1-iysub)/2.0))-1.0)*5 +1)*nyp/36.0+0.5));
1201 0 : if(linear)
1202 0 : nysub=Int(floor((ceil(fabs(float(2*elrow+1-iysub)/2.0))/factor)*nyp + 0.5));
1203 : else
1204 0 : nysub=Int(floor((sqrt(ceil(fabs(float(2*elrow+1-iysub)/2.0)))/factor)*nyp + 0.5));
1205 : //nysub=Int(floor((ceil(fabs(float(2*elrow+1-iysub)/2.0))*ceil(fabs(float(2*elrow+1-iysub)/2.0))*ceil(fabs(float(2*elrow+1-iysub)/2.0))/factor)*nyp + 0.5));
1206 0 : y0=miny;
1207 0 : elrow-=1;
1208 :
1209 0 : while(elrow >=0){
1210 : //y0+=Int(floor(((ceil(fabs(float(2*elrow+1-iysub)/2.0))-1.0)*5 +1)*nyp/36.0+0.5));
1211 0 : if(linear)
1212 0 : y0+=Int(floor((ceil(fabs(float(2*elrow+1-iysub)/2.0))/factor)*nyp + 0.5));
1213 : else
1214 0 : y0+=Int(floor((sqrt(ceil(fabs(float(2*elrow+1-iysub)/2.0)))/factor)*nyp + 0.5));
1215 : //y0+=Int(floor((ceil(fabs(float(2*elrow+1-iysub)/2.0))*ceil(fabs(float(2*elrow+1-iysub)/2.0))*ceil(fabs(float(2*elrow+1-iysub)/2.0))/factor)*nyp + 0.5));
1216 0 : elrow-=1;
1217 : }
1218 : }
1219 :
1220 :
1221 0 : y0+=1;
1222 0 : x0+=1;
1223 :
1224 :
1225 0 : }
1226 :
1227 0 : void WProjectFT::tweakGridSector(const Int& nx, const Int& ny, const Int& ixsub, const Int& iysub, Vector<Int>& xsect, Vector<Int>& ysect, Vector<Int>& nxsect, Vector<Int>& nysect){
1228 0 : Vector<Int> x0, y0, nxsub, nysub;
1229 0 : Vector<Float> xcut(nx/2);
1230 0 : Vector<Float> ycut(ny/2);
1231 0 : if(griddedData2.nelements() >0 ){
1232 : //cerr << "shapes " << xcut.shape() << " gd " << amplitude(griddedData2(IPosition(4, 0, ny/2-1, 0, 0), IPosition(4, nx/2-1, ny/2-1, 0,0))).shape() << endl;
1233 0 : convertArray(xcut, amplitude(griddedData2(IPosition(4, 0, ny/2-1, 0, 0), IPosition(4, nx/2-1, ny/2-1, 0,0))).reform(xcut.shape()));
1234 0 : convertArray(ycut, amplitude(griddedData2(IPosition(4, nx/2-1, 0, 0, 0), IPosition(4, nx/2-1, ny/2-1, 0,0))).reform(ycut.shape()));
1235 : }
1236 : else{
1237 0 : xcut=amplitude(griddedData(IPosition(4, 0, ny/2-1, 0, 0), IPosition(4, nx/2-1, ny/2-1, 0,0)));
1238 0 : ycut=amplitude(griddedData(IPosition(4, nx/2-1, 0, 0, 0), IPosition(4, nx/2-1, ny/2-1, 0,0)));
1239 : }
1240 : //cerr << griddedData2.shape() << " " << griddedData.shape() << endl;
1241 0 : Vector<Float> cumSumX(nx/2, 0);
1242 : //Vector<Float> cumSumX2(nx/2,0);
1243 0 : cumSumX(0)=xcut(0);
1244 : //cumSumX2(0)=cumSumX(0)*cumSumX(0);
1245 0 : Float sumX=sum(xcut);
1246 0 : if(sumX==0.0)
1247 0 : return;
1248 : //cerr << "sumX " << sumX << endl;
1249 : //sumX *= sumX;
1250 0 : x0.resize(ixsub);
1251 0 : x0=nx/2-1;
1252 0 : nxsub.resize(ixsub);
1253 0 : nxsub=0;
1254 0 : x0(0)=0;
1255 0 : Int counter=1;
1256 0 : for (Int k=1; k < nx/2; ++k){
1257 0 : cumSumX(k)=cumSumX(k-1)+xcut(k);
1258 : //cumSumX2(k)=cumSumX(k)*cumSumX(k);
1259 0 : Float nextEdge=sumX/(Float(ixsub/2)*Float(ixsub/2))*Float(counter)*Float(counter);
1260 0 : if(cumSumX(k-1) < nextEdge && cumSumX(k) >= nextEdge){
1261 0 : x0(counter)=k;
1262 : //cerr << counter << " " << k << " diff " << x0(counter)-x0(counter-1) << endl;
1263 0 : nxsub(counter-1)=x0(counter)-x0(counter-1);
1264 0 : ++counter;
1265 : }
1266 : }
1267 :
1268 0 : x0(ixsub/2)=nx/2-1;
1269 0 : nxsub(ixsub/2)=nxsub(ixsub/2-1)=nx/2-1-x0(ixsub/2-1);
1270 0 : for(Int k=ixsub/2+1; k < ixsub; ++k){
1271 0 : x0(k)=x0(k-1)+ nxsub(ixsub-k);
1272 0 : nxsub(k)=nxsub(ixsub-k-1);
1273 : }
1274 0 : nxsub(ixsub-1)+=1;
1275 :
1276 0 : Vector<Float> cumSumY(ny/2, 0);
1277 : //Vector<Float> cumSumY2(ny/2,0);
1278 0 : cumSumY(0)=ycut(0);
1279 : //cumSumY2(0)=cumSumY(0)*cumSumY(0);
1280 0 : Float sumY=sum(ycut);
1281 0 : if(sumY==0.0)
1282 0 : return;
1283 : //sumY *=sumY;
1284 0 : y0.resize(iysub);
1285 0 : y0=ny/2-1;
1286 0 : nysub.resize(iysub);
1287 0 : nysub=0;
1288 0 : y0(0)=0;
1289 0 : counter=1;
1290 0 : for (Int k=1; k < ny/2; ++k){
1291 0 : cumSumY(k)=cumSumY(k-1)+ycut(k);
1292 : //cumSumY2(k)=cumSumY(k)*cumSumY(k);
1293 0 : Float nextEdge=sumY/(Float(iysub/2)*Float(iysub/2))*Float(counter)*Float(counter);
1294 0 : if(cumSumY(k-1) < nextEdge && cumSumY(k) >= nextEdge){
1295 0 : y0(counter)=k;
1296 0 : nysub(counter-1)=y0(counter)-y0(counter-1);
1297 0 : ++counter;
1298 : }
1299 : }
1300 :
1301 0 : y0(ixsub/2)=ny/2-1;
1302 0 : nysub(iysub/2)=nysub(iysub/2-1)=ny/2-1-y0(iysub/2-1);
1303 0 : for(Int k=iysub/2+1; k < iysub; ++k){
1304 0 : y0(k)=y0(k-1)+ nysub(iysub-k);
1305 0 : nysub(k)=nysub(iysub-k-1);
1306 : }
1307 0 : nysub(iysub-1)+=1;
1308 :
1309 :
1310 : //cerr << " x0 " << x0 << " nxsub " << nxsub << endl;
1311 : //cerr << " y0 " << y0 << " nysub " << nysub << endl;
1312 0 : x0+=1;
1313 0 : y0+=1;
1314 0 : xsect.resize(ixsub*iysub);
1315 0 : ysect.resize(ixsub*iysub);
1316 0 : nxsect.resize(ixsub*iysub);
1317 0 : nysect.resize(ixsub*iysub);
1318 0 : for (Int iy=0; iy < iysub; ++iy){
1319 0 : for (Int ix=0; ix< ixsub; ++ix){
1320 :
1321 0 : xsect(iy*ixsub+ix)=x0[ix];
1322 0 : ysect(iy*ixsub+ix)=y0[iy];
1323 0 : nxsect(iy*ixsub+ix)=nxsub[ix];
1324 0 : nysect(iy*ixsub+ix)=nysub[iy];
1325 : }
1326 : }
1327 :
1328 :
1329 :
1330 : }
1331 :
1332 :
1333 0 : void WProjectFT::get(VisBuffer& vb, Int row)
1334 : {
1335 :
1336 : ///This is needed here as virtual model column may not call initializeToVis
1337 0 : findConvFunction(*image, vb);
1338 : // If row is -1 then we pass through all rows
1339 : Int startRow, endRow, nRow;
1340 0 : if (row==-1) {
1341 0 : nRow=vb.nRow();
1342 0 : startRow=0;
1343 0 : endRow=nRow-1;
1344 : //vb.modelVisCube()=Complex(0.0,0.0);
1345 : } else {
1346 0 : nRow=1;
1347 0 : startRow=row;
1348 0 : endRow=row;
1349 : //vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
1350 : }
1351 :
1352 : // Get the uvws in a form that Fortran can use
1353 0 : Matrix<Double> uvw(3, vb.uvw().nelements());
1354 0 : uvw=0.0;
1355 0 : Vector<Double> dphase(vb.uvw().nelements());
1356 0 : dphase=0.0;
1357 : //NEGATING to correct for an image inversion problem
1358 0 : for (Int i=startRow;i<=endRow;i++) {
1359 0 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
1360 0 : uvw(2,i)=vb.uvw()(i)(2);
1361 : }
1362 :
1363 0 : rotateUVW(uvw, dphase, vb);
1364 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1365 :
1366 : // This is the convention for dphase
1367 : // dphase*=-1.0;
1368 :
1369 :
1370 : //Check if ms has changed then cache new spw and chan selection
1371 0 : if(vb.newMS())
1372 0 : matchAllSpwChans(vb);
1373 : //Here we redo the match or use previous match
1374 :
1375 : //Channel matching for the actual spectral window of buffer
1376 0 : if(doConversion_p[vb.spectralWindow()]){
1377 0 : matchChannel(vb.spectralWindow(), vb);
1378 : }
1379 : else{
1380 0 : chanMap.resize();
1381 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
1382 : }
1383 :
1384 : //No point in reading data if its not matching in frequency
1385 0 : if(max(chanMap)==-1)
1386 0 : return;
1387 :
1388 0 : Cube<Complex> data;
1389 0 : Cube<Int> flags;
1390 0 : getInterpolateArrays(vb, data, flags);
1391 :
1392 :
1393 :
1394 : Complex *datStorage;
1395 : Bool isCopy;
1396 0 : datStorage=data.getStorage(isCopy);
1397 :
1398 0 : Vector<Int> rowFlags(vb.nRow());
1399 0 : rowFlags=0;
1400 0 : rowFlags(vb.flagRow())=true;
1401 0 : if(!usezero_p) {
1402 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1403 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
1404 : }
1405 : }
1406 :
1407 0 : Int nvp=data.shape()(0);
1408 0 : Int nvc=data.shape()(1);
1409 0 : Int nvisrow=data.shape()(2);
1410 0 : Int nc=nchan;
1411 0 : Int np=npol;
1412 0 : Int nxp=nx;
1413 0 : Int nyp=ny;
1414 0 : Cube<Int> loc(3, nvc, nvisrow);
1415 0 : Cube<Int> off(3, nvc, nRow);
1416 0 : Int csamp=convSampling;
1417 0 : Int csize=convSize;
1418 0 : Int wcsize=wConvSize;
1419 0 : Matrix<Complex> phasor(nvc, nRow);
1420 : Bool delphase;
1421 : Bool del;
1422 0 : Complex * phasorstor=phasor.getStorage(delphase);
1423 0 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1424 0 : const Double * scalestor=uvScale.getStorage(del);
1425 0 : const Double * offsetstor=uvOffset.getStorage(del);
1426 0 : Int * locstor=loc.getStorage(del);
1427 0 : Int * offstor=off.getStorage(del);
1428 0 : const Int * flagstor=flags.getStorage(del);
1429 0 : const Int * rowflagstor=rowFlags.getStorage(del);
1430 0 : const Double *dpstor=dphase.getStorage(del);
1431 : Bool uvwcopy;
1432 0 : const Double *uvwstor=uvw.getStorage(uvwcopy);
1433 : Bool gridcopy;
1434 0 : const Complex *gridstor=griddedData.getStorage(gridcopy);
1435 : Bool convcopy;
1436 0 : const Complex *convstor=convFunc.getStorage(convcopy);
1437 0 : const Int* pmapstor=polMap.getStorage(del);
1438 0 : const Int *cmapstor=chanMap.getStorage(del);
1439 0 : const Int * suppstor=convSupport.getStorage(del);
1440 : Int irow;
1441 0 : Int nth=1;
1442 : #ifdef _OPENMP
1443 0 : if(numthreads_p >0){
1444 0 : nth=min(numthreads_p, omp_get_max_threads());
1445 : }
1446 : else{
1447 0 : nth= omp_get_max_threads();
1448 : }
1449 : // nth=min(4,nth);
1450 : #endif
1451 0 : Int dow=1;
1452 0 : Double cinv=Double(1.0)/C::c;
1453 :
1454 0 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvc, scalestor, offsetstor, csamp, phasorstor, uvwstor, locstor, offstor, dpstor, dow, cinv) shared(startRow, endRow) num_threads(nth)
1455 : {
1456 : #pragma omp for schedule(dynamic)
1457 : for (irow=startRow; irow<=endRow; ++irow){
1458 : /*locateuvw(uvwstor,dpstor, visfreqstor, nvc, scalestor, offsetstor, csamp,
1459 : locstor,
1460 : offstor, phasorstor, irow, true);*/
1461 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
1462 : }
1463 :
1464 : }//end pragma parallel
1465 0 : Int rbeg=startRow+1;
1466 0 : Int rend=endRow+1;
1467 0 : Int npart=nth*2;
1468 :
1469 0 : Timer tim;
1470 0 : tim.mark();
1471 :
1472 0 : Int ix=0;
1473 0 : #pragma omp parallel default(none) private(ix, rbeg, rend) firstprivate(uvwstor, datStorage, flagstor, rowflagstor, convstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc, suppstor, csamp, csize, wcsize, nvp, nvc, nvisrow, phasorstor, locstor, offstor) shared(npart) num_threads(nth)
1474 : {
1475 : #pragma omp for schedule(dynamic,1)
1476 : for (ix=0; ix< npart; ++ix){
1477 : rbeg=ix*(nvisrow/npart)+1;
1478 : rend=(ix != (npart-1)) ? (rbeg+(nvisrow/npart)-1) : (rbeg+(nvisrow/npart)+nvisrow%npart-1) ;
1479 : // cerr << "rbeg " << rbeg << " rend " << rend << " nRow " << nvisrow << endl;
1480 : sectdwgrid(uvwstor,
1481 : datStorage,
1482 : &nvp,
1483 : &nvc,
1484 : flagstor,
1485 : rowflagstor,
1486 : &nvisrow,
1487 : gridstor,
1488 : &nxp,
1489 : &nyp,
1490 : &np,
1491 : &nc,
1492 : suppstor,
1493 : &csize,
1494 : &csamp,
1495 : &wcsize,
1496 : convstor,
1497 : cmapstor,
1498 : pmapstor,
1499 : &rbeg, &rend, locstor, offstor, phasorstor);
1500 : }
1501 :
1502 : }//end pragma parallel
1503 0 : data.putStorage(datStorage, isCopy);
1504 0 : uvw.freeStorage(uvwstor, uvwcopy);
1505 0 : griddedData.freeStorage(gridstor, gridcopy);
1506 0 : convFunc.freeStorage(convstor, convcopy);
1507 0 : timedegrid_p+=tim.real();
1508 :
1509 0 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1510 : }
1511 :
1512 :
1513 :
1514 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1515 : // return the resulting image
1516 0 : ImageInterface<Complex>& WProjectFT::getImage(Matrix<Float>& weights,
1517 : Bool normalize)
1518 : {
1519 : //AlwaysAssert(lattice, AipsError);
1520 0 : AlwaysAssert(image, AipsError);
1521 :
1522 0 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
1523 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
1524 : }
1525 0 : logIO() << LogOrigin("WProjectFT", "getImage") << LogIO::NORMAL;
1526 :
1527 0 : weights.resize(sumWeight.shape());
1528 :
1529 0 : convertArray(weights, sumWeight);
1530 :
1531 : // If the weights are all zero then we cannot normalize
1532 : // otherwise we don't care.
1533 0 : if(max(weights)==0.0) {
1534 0 : if(normalize) {
1535 : logIO() << LogIO::SEVERE
1536 : << "No useful data in WProjectFT: weights all zero"
1537 0 : << LogIO::POST;
1538 : }
1539 : else {
1540 : logIO() << LogIO::WARN
1541 : << "No useful data in WProjectFT: weights all zero"
1542 0 : << LogIO::POST;
1543 : }
1544 : }
1545 : else {
1546 : logIO() << LogIO::DEBUGGING
1547 0 : << "Starting FFT and scaling of image" << LogIO::POST;
1548 :
1549 0 : if(useDoubleGrid_p){
1550 0 : ArrayLattice<DComplex> darrayLattice(griddedData2);
1551 0 : LatticeFFT::cfft2d(darrayLattice,false);
1552 0 : griddedData.resize(griddedData2.shape());
1553 0 : convertArray(griddedData, griddedData2);
1554 0 : SynthesisUtilMethods::getResource("mem peak in getImage");
1555 0 : griddedData2.resize();
1556 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1557 0 : lattice=arrayLattice;
1558 : }else{
1559 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1560 0 : lattice=arrayLattice;
1561 0 : LatticeFFT::cfft2d(*lattice,false);
1562 :
1563 : }
1564 :
1565 :
1566 0 : const IPosition latticeShape = lattice->shape();
1567 :
1568 :
1569 : {
1570 0 : Int inx = lattice->shape()(0);
1571 0 : Int iny = lattice->shape()(1);
1572 0 : Vector<Complex> correction(inx);
1573 0 : correction=Complex(1.0, 0.0);
1574 0 : Vector<Float> sincConvX(nx);
1575 0 : for (Int ix=0;ix<nx;ix++) {
1576 0 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
1577 0 : if(ix==nx/2) {
1578 0 : sincConvX(ix)=1.0;
1579 : }
1580 : else {
1581 0 : sincConvX(ix)=sin(x)/x;
1582 : }
1583 : }
1584 0 : Vector<Float> sincConvY(ny);
1585 0 : for (Int ix=0;ix<ny;ix++) {
1586 0 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
1587 0 : if(ix==ny/2) {
1588 0 : sincConvY(ix)=1.0;
1589 : }
1590 : else {
1591 0 : sincConvY(ix)=sin(x)/x;
1592 : }
1593 : }
1594 :
1595 :
1596 :
1597 0 : IPosition cursorShape(4, inx, 1, 1, 1);
1598 0 : IPosition axisPath(4, 0, 1, 2, 3);
1599 0 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1600 0 : LatticeIterator<Complex> lix(*lattice, lsx);
1601 0 : for(lix.reset();!lix.atEnd();lix++) {
1602 0 : Int pol=lix.position()(2);
1603 0 : Int chan=lix.position()(3);
1604 0 : if(weights(pol, chan)!=0.0) {
1605 0 : Int iy=lix.position()(1);
1606 0 : gridder->correctX1D(correction,iy);
1607 0 : for (Int ix=0;ix<nx;ix++) {
1608 0 : correction(ix)*=sincConvX(ix)*sincConvY(iy);
1609 : }
1610 0 : lix.rwVectorCursor()/=correction;
1611 0 : if(normalize) {
1612 0 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
1613 0 : lix.rwCursor()*=rnorm;
1614 : }
1615 : else {
1616 0 : Complex rnorm(Float(inx)*Float(iny));
1617 0 : lix.rwCursor()*=rnorm;
1618 : }
1619 : }
1620 : else {
1621 0 : lix.woCursor()=0.0;
1622 : }
1623 : }
1624 : }
1625 :
1626 0 : if(!isTiled) {
1627 : // Check the section from the image BEFORE converting to a lattice
1628 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
1629 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1630 0 : IPosition stride(4, 1);
1631 0 : IPosition trc(blc+image->shape()-stride);
1632 :
1633 : // Do the copy
1634 0 : image->put(griddedData(blc, trc));
1635 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
1636 0 : griddedData.resize(IPosition(1,0));
1637 : }
1638 : }
1639 0 : if(!arrayLattice.null()) arrayLattice=0;
1640 0 : if(!lattice.null()) lattice=0;
1641 0 : griddedData.resize();
1642 0 : return *image;
1643 : }
1644 :
1645 : // Get weight image
1646 0 : void WProjectFT::getWeightImage(ImageInterface<Float>& weightImage,
1647 : Matrix<Float>& weights)
1648 : {
1649 :
1650 0 : logIO() << LogOrigin("WProjectFT", "getWeightImage") << LogIO::NORMAL;
1651 :
1652 0 : weights.resize(sumWeight.shape());
1653 0 : convertArray(weights,sumWeight);
1654 :
1655 0 : const IPosition latticeShape = weightImage.shape();
1656 :
1657 0 : Int inx=latticeShape(0);
1658 0 : Int iny=latticeShape(1);
1659 :
1660 0 : IPosition loc(2, 0);
1661 0 : IPosition cursorShape(4, inx, iny, 1, 1);
1662 0 : IPosition axisPath(4, 0, 1, 2, 3);
1663 0 : LatticeStepper lsx(latticeShape, cursorShape, axisPath);
1664 0 : LatticeIterator<Float> lix(weightImage, lsx);
1665 0 : for(lix.reset();!lix.atEnd();lix++) {
1666 0 : Int pol=lix.position()(2);
1667 0 : Int chan=lix.position()(3);
1668 0 : lix.rwCursor()=weights(pol,chan);
1669 : }
1670 0 : }
1671 :
1672 0 : Bool WProjectFT::toRecord(String& error,
1673 : RecordInterface& outRec, Bool withImage, const String diskimage)
1674 : {
1675 :
1676 : /*
1677 :
1678 : CountedPtr<WPConvFunc> wpConvFunc_p;
1679 : */
1680 :
1681 : // Save the current WProjectFT object to an output state record
1682 0 : Bool retval = true;
1683 : //save the base class variables
1684 : //this is a memory hog and slow on saving and recovering...better to recompute convfunctions
1685 : /* Record wpconvrec;
1686 : if(wpConvFunc_p->toRecord(wpconvrec))
1687 : outRec.defineRecord("wpconvfunc", wpconvrec);
1688 : */
1689 0 : Float elpadd=padding_p;
1690 0 : if(toVis_p && withImage)
1691 0 : elpadd=1.0;
1692 0 : if(!FTMachine::toRecord(error, outRec, withImage, diskimage))
1693 0 : return false;
1694 :
1695 0 : outRec.define("uvscale", uvScale);
1696 0 : outRec.define("uvoffset", uvOffset);
1697 0 : outRec.define("istiled", isTiled);
1698 0 : outRec.define("cachesize", Int64(cachesize));
1699 0 : outRec.define("tilesize", tilesize);
1700 0 : outRec.define("maxabsdata", maxAbsData);
1701 0 : Vector<Int> center_loc(4), offset_loc(4);
1702 0 : for (Int k=0; k<4 ; k++){
1703 0 : center_loc(k)=centerLoc(k);
1704 0 : offset_loc(k)=offsetLoc(k);
1705 : }
1706 0 : outRec.define("centerloc", center_loc);
1707 0 : outRec.define("offsetloc", offset_loc);
1708 0 : outRec.define("padding", elpadd);
1709 0 : outRec.define("nwplanes", nWPlanes_p);
1710 0 : outRec.define("savedwscale", savedWScale_p);
1711 0 : outRec.define("usezero", usezero_p);
1712 : ///no need to save convfunc as it can be big and is recalculated anyways
1713 : ///outRec.define("convfunc", convFunc);
1714 0 : outRec.define("convsampling", convSampling);
1715 0 : outRec.define("convsize", convSize);
1716 0 : outRec.define("convsupport", convSupport);
1717 0 : outRec.define("convsizes", convSizes_p);
1718 0 : outRec.define("wconvsize", wConvSize);
1719 0 : outRec.define("lastindex", lastIndex_p);
1720 0 : outRec.define("minw", minW_p);
1721 0 : outRec.define("maxw", maxW_p);
1722 0 : outRec.define("rmsw", rmsW_p);
1723 :
1724 :
1725 0 : return retval;
1726 : }
1727 :
1728 0 : Bool WProjectFT::fromRecord(String& error,
1729 : const RecordInterface& inRec)
1730 : {
1731 0 : if(!FTMachine::fromRecord(error, inRec))
1732 0 : return false;
1733 0 : machineName_p="WProjectFT";
1734 0 : Bool retval = true;
1735 0 : imageCache=0; lattice=0; arrayLattice=0;
1736 0 : inRec.get("uvscale", uvScale);
1737 0 : inRec.get("uvoffset", uvOffset);
1738 0 : inRec.get("istiled", isTiled);
1739 0 : cachesize=inRec.asInt64("cachesize");
1740 0 : inRec.get("tilesize", tilesize);
1741 0 : inRec.get("maxabsdata", maxAbsData);
1742 0 : Vector<Int> center_loc(4), offset_loc(4);
1743 0 : inRec.get("centerloc", center_loc);
1744 0 : inRec.get("offsetloc", offset_loc);
1745 0 : uInt ndim4 = 4;
1746 0 : centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2),
1747 0 : center_loc(3));
1748 0 : offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2),
1749 0 : offset_loc(3));
1750 0 : inRec.get("minw", minW_p);
1751 0 : inRec.get("maxw", maxW_p);
1752 0 : inRec.get("rmsw", rmsW_p);
1753 0 : if(inRec.isDefined("wpconvfunc")){
1754 0 : wpConvFunc_p=new WPConvFunc(inRec.asRecord("wpconvfunc"));
1755 : }
1756 : else{
1757 0 : wpConvFunc_p=new WPConvFunc(minW_p, maxW_p, rmsW_p);
1758 : }
1759 0 : inRec.get("padding", padding_p);
1760 0 : inRec.get("nwplanes", nWPlanes_p);
1761 0 : inRec.get("savedwscale", savedWScale_p);
1762 0 : inRec.get("usezero", usezero_p);
1763 : ///have stopped saving this parameter
1764 : ////inRec.get("convfunc", convFunc);
1765 0 : convFunc.resize();
1766 0 : inRec.get("convsampling", convSampling);
1767 0 : inRec.get("convsize", convSize);
1768 0 : inRec.get("convsupport", convSupport);
1769 0 : inRec.get("convsizes", convSizes_p);
1770 0 : inRec.get("wconvsize", wConvSize);
1771 0 : inRec.get("lastindex", lastIndex_p);
1772 :
1773 0 : gridder=0;
1774 : ///setup some of the parameters
1775 0 : init();
1776 :
1777 :
1778 :
1779 0 : return retval;
1780 : }
1781 :
1782 0 : void WProjectFT::ok() {
1783 0 : AlwaysAssert(image, AipsError);
1784 0 : }
1785 :
1786 : // Make a plain straightforward honest-to-God image. This returns
1787 : // a complex image, without conversion to Stokes. The representation
1788 : // is that required for the visibilities.
1789 : //----------------------------------------------------------------------
1790 0 : void WProjectFT::makeImage(FTMachine::Type type,
1791 : VisSet& vs,
1792 : ImageInterface<Complex>& theImage,
1793 : Matrix<Float>& weight) {
1794 :
1795 :
1796 0 : logIO() << LogOrigin("WProjectFT", "makeImage") << LogIO::NORMAL;
1797 :
1798 0 : if(type==FTMachine::COVERAGE) {
1799 : logIO() << "Type COVERAGE not defined for Fourier transforms"
1800 0 : << LogIO::EXCEPTION;
1801 : }
1802 :
1803 :
1804 : // Initialize the gradients
1805 0 : ROVisIter& vi(vs.iter());
1806 :
1807 : // Loop over all visibilities and pixels
1808 0 : VisBuffer vb(vi);
1809 :
1810 : // Initialize put (i.e. transform to Sky) for this model
1811 0 : vi.origin();
1812 :
1813 0 : if(vb.polFrame()==MSIter::Linear) {
1814 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
1815 : }
1816 : else {
1817 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
1818 : }
1819 :
1820 0 : initializeToSky(theImage,weight,vb);
1821 :
1822 : // Loop over the visibilities, putting VisBuffers
1823 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1824 0 : for (vi.origin(); vi.more(); vi++) {
1825 :
1826 0 : switch(type) {
1827 0 : case FTMachine::RESIDUAL:
1828 0 : vb.visCube()=vb.correctedVisCube();
1829 0 : vb.visCube()-=vb.modelVisCube();
1830 0 : put(vb, -1, false);
1831 0 : break;
1832 0 : case FTMachine::MODEL:
1833 0 : vb.visCube()=vb.modelVisCube();
1834 0 : put(vb, -1, false);
1835 0 : break;
1836 0 : case FTMachine::CORRECTED:
1837 0 : vb.visCube()=vb.correctedVisCube();
1838 0 : put(vb, -1, false);
1839 0 : break;
1840 0 : case FTMachine::PSF:
1841 0 : vb.visCube()=Complex(1.0,0.0);
1842 0 : put(vb, -1, true);
1843 0 : break;
1844 0 : case FTMachine::OBSERVED:
1845 : default:
1846 0 : put(vb, -1, false);
1847 0 : break;
1848 : }
1849 : }
1850 : }
1851 0 : finalizeToSky();
1852 : // Normalize by dividing out weights, etc.
1853 0 : getImage(weight, true);
1854 0 : }
1855 :
1856 :
1857 0 : String WProjectFT::name() const {
1858 :
1859 0 : return machineName_p;
1860 :
1861 : }
1862 :
1863 0 : void WProjectFT::wStat(ROVisibilityIterator& vi, Double& minW, Double& maxW, Double& rmsW){
1864 0 : VisBuffer vb(vi);
1865 0 : maxW=0.0;
1866 0 : minW=1e99;
1867 0 : Double nval=0;
1868 0 : rmsW=0.0;
1869 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
1870 0 : for (vi.origin();vi.more();vi++) {
1871 0 : maxW=max(maxW, max(abs(vb.uvwMat().row(2)*max(vb.frequency())))/C::c);
1872 0 : minW=min(minW, min(abs(vb.uvwMat().row(2)*min(vb.frequency())))/C::c);
1873 : ///////////some shenanigans as some compilers is confusing * operator for vector
1874 0 : Vector<Double> elw;
1875 0 : elw=vb.uvwMat().row(2);
1876 0 : elw*=vb.uvwMat().row(2);
1877 : //////////////////////////////////////////////////
1878 0 : rmsW+=sum(elw);
1879 :
1880 0 : nval+=Double(vb.nRow());
1881 : }
1882 : }
1883 0 : rmsW=sqrt(rmsW/Double(nval))*max(vb.frequency())/C::c;
1884 :
1885 0 : }
1886 :
1887 :
1888 : } //# NAMESPACE CASA - END
1889 :
|
