Line data Source code
1 : //# MosaicFT.cc: Implementation of MosaicFT class
2 : //# Copyright (C) 2002-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 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: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Quanta/UnitMap.h>
29 : #include <casacore/casa/Quanta/MVTime.h>
30 : #include <casacore/casa/Quanta/UnitVal.h>
31 : #include <casacore/measures/Measures/Stokes.h>
32 : #include <casacore/measures/Measures/UVWMachine.h>
33 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
34 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
35 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
36 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
37 : #include <casacore/coordinates/Coordinates/Projection.h>
38 : #include <casacore/ms/MeasurementSets/MSColumns.h>
39 : #include <casacore/casa/BasicSL/Constants.h>
40 : #include <casacore/scimath/Mathematics/FFTServer.h>
41 : #include <synthesis/TransformMachines2/MosaicFT.h>
42 : #include <synthesis/TransformMachines2/SimplePBConvFunc.h>
43 : #include <synthesis/TransformMachines2/HetArrayConvFunc.h>
44 : #include <synthesis/TransformMachines/PBMath.h>
45 : #include <synthesis/TransformMachines2/VPSkyJones.h>
46 : #include <casacore/scimath/Mathematics/RigidVector.h>
47 : #include <msvis/MSVis/StokesVector.h>
48 : #include <synthesis/TransformMachines/StokesImageUtil.h>
49 : #include <msvis/MSVis/VisBuffer2.h>
50 : #include <msvis/MSVis/VisibilityIterator2.h>
51 : #include <casacore/images/Images/ImageInterface.h>
52 : #include <casacore/images/Images/PagedImage.h>
53 : #include <casacore/images/Images/SubImage.h>
54 : #include <casacore/images/Regions/ImageRegion.h>
55 : #include <casacore/images/Regions/WCBox.h>
56 : #include <casacore/casa/Containers/Block.h>
57 : #include <casacore/casa/Containers/Record.h>
58 : #include <casacore/casa/Arrays/ArrayLogical.h>
59 : #include <casacore/casa/Arrays/ArrayMath.h>
60 : #include <casacore/casa/Arrays/Array.h>
61 : #include <casacore/casa/Arrays/MaskedArray.h>
62 : #include <casacore/casa/Arrays/Vector.h>
63 : #include <casacore/casa/Arrays/Slice.h>
64 : #include <casacore/casa/Arrays/Matrix.h>
65 : #include <casacore/casa/Arrays/Cube.h>
66 : #include <casacore/casa/Arrays/MatrixIter.h>
67 : #include <casacore/casa/BasicSL/String.h>
68 : #include <casacore/casa/Utilities/Assert.h>
69 : #include <casacore/casa/Exceptions/Error.h>
70 : #include <casacore/lattices/Lattices/ArrayLattice.h>
71 : #include <casacore/lattices/Lattices/SubLattice.h>
72 : #include <casacore/lattices/LRegions/LCBox.h>
73 : #include <casacore/lattices/LEL/LatticeExpr.h>
74 : #include <casacore/lattices/Lattices/LatticeCache.h>
75 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
76 : #include <casacore/lattices/Lattices/LatticeIterator.h>
77 : #include <casacore/lattices/Lattices/LatticeStepper.h>
78 : #include <casacore/casa/Utilities/CompositeNumber.h>
79 : #include <casacore/casa/OS/Timer.h>
80 : #include <casacore/casa/OS/HostInfo.h>
81 : #include <sstream>
82 : #ifdef _OPENMP
83 : #include <omp.h>
84 : #endif
85 : using namespace casacore;
86 : namespace casa { //# NAMESPACE CASA - BEGIN
87 : namespace refim {//# namespace for imaging refactor
88 : using namespace casacore;
89 : using namespace casa;
90 : using namespace casacore;
91 : using namespace casa::refim;
92 :
93 334 : MosaicFT::MosaicFT(SkyJones* sj, MPosition mloc, String stokes,
94 : Long icachesize, Int itilesize,
95 334 : Bool usezero, Bool useDoublePrec, Bool useConjConvFunc, Bool usePointing)
96 334 : : FTMachine(), sj_p(sj),
97 334 : imageCache(0), cachesize(icachesize), tilesize(itilesize), gridder(nullptr),
98 334 : isTiled(false),
99 334 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
100 1336 : mspc(0), msac(0), pointingToImage(0), usezero_p(usezero), convSampling(1),
101 1002 : skyCoverage_p( ), machineName_p("MosaicFT"), stokes_p(stokes), useConjConvFunc_p(useConjConvFunc), usePointingTable_p(usePointing),timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0)
102 : {
103 334 : convSize=0;
104 334 : lastIndex_p=0;
105 334 : doneWeightImage_p=false;
106 334 : convWeightImage_p=0;
107 334 : pbConvFunc_p=new SimplePBConvFunc();
108 :
109 334 : mLocation_p=mloc;
110 334 : useDoubleGrid_p=useDoublePrec;
111 : // We should get rid of the ms dependence in the constructor
112 : // not used
113 334 : }
114 :
115 0 : MosaicFT::MosaicFT(const RecordInterface& stateRec)
116 0 : : FTMachine()
117 : {
118 : // Construct from the input state record
119 0 : String error;
120 0 : if (!fromRecord(error, stateRec)) {
121 0 : throw (AipsError("Failed to create MosaicFT: " + error));
122 : };
123 0 : }
124 :
125 : //----------------------------------------------------------------------
126 297 : MosaicFT& MosaicFT::operator=(const MosaicFT& other)
127 : {
128 297 : if(this!=&other) {
129 :
130 : //Do the base parameters
131 297 : FTMachine::operator=(other);
132 :
133 297 : convSampling=other.convSampling;
134 297 : sj_p=other.sj_p;
135 297 : imageCache=other.imageCache;
136 297 : cachesize=other.cachesize;
137 297 : tilesize=other.tilesize;
138 297 : isTiled=other.isTiled;
139 : //lattice=other.lattice;
140 297 : lattice=0;
141 : // arrayLattice=other.arrayLattice;
142 : // weightLattice=other.weightLattice;
143 : //if(arrayLattice) delete arrayLattice;
144 297 : arrayLattice=0;
145 : //if(weightLattice) delete weightLattice;
146 297 : weightLattice=0;
147 297 : maxAbsData=other.maxAbsData;
148 297 : centerLoc=other.centerLoc;
149 297 : offsetLoc=other.offsetLoc;
150 297 : pointingToImage=other.pointingToImage;
151 297 : usezero_p=other.usezero_p;
152 297 : doneWeightImage_p=other.doneWeightImage_p;
153 297 : pbConvFunc_p=other.pbConvFunc_p;
154 297 : stokes_p=other.stokes_p;
155 297 : if(!other.skyCoverage_p.null())
156 0 : skyCoverage_p=other.skyCoverage_p;
157 : else
158 297 : skyCoverage_p=0;
159 297 : if(other.convWeightImage_p !=0)
160 0 : convWeightImage_p=(TempImage<Complex> *)other.convWeightImage_p->cloneII();
161 : else
162 297 : convWeightImage_p=0;
163 297 : if(other.gridder==0)
164 297 : gridder.reset(nullptr);
165 : else{
166 0 : uvScale=other.uvScale;
167 0 : uvOffset=other.uvOffset;
168 0 : gridder.reset(new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
169 0 : uvScale, uvOffset,
170 0 : "SF"));
171 :
172 : }
173 297 : useConjConvFunc_p=other.useConjConvFunc_p;
174 297 : usePointingTable_p=other.usePointingTable_p;
175 297 : timemass_p=other.timemass_p;
176 297 : timegrid_p=other.timegrid_p;
177 297 : timedegrid_p=other.timedegrid_p;
178 : };
179 297 : return *this;
180 : };
181 :
182 : //----------------------------------------------------------------------
183 294 : MosaicFT::MosaicFT(const MosaicFT& other): FTMachine(),machineName_p("MosaicFT")
184 : {
185 294 : operator=(other);
186 294 : }
187 :
188 : //----------------------------------------------------------------------
189 : //void MosaicFT::setSharingFT(MosaicFT& otherFT){
190 : // otherFT_p=&otherFT;
191 : //}
192 491 : void MosaicFT::init(const vi::VisBuffer2& /*vb*/) {
193 :
194 : /* if((image->shape().product())>cachesize) {
195 : isTiled=true;
196 : }
197 : else {
198 : isTiled=false;
199 : }
200 : */
201 : //For now only isTiled false works.
202 491 : isTiled=false;
203 491 : nx = image->shape()(0);
204 491 : ny = image->shape()(1);
205 491 : npol = image->shape()(2);
206 491 : nchan = image->shape()(3);
207 :
208 :
209 : // if(skyCoverage_p==0){
210 : // Double memoryMB=HostInfo::memoryTotal()/1024.0/(20.0);
211 : // skyCoverage_p=new TempImage<Float> (IPosition(4,nx,ny,1,1),
212 : // image->coordinates(), memoryMB);
213 : //Setting it to zero
214 : // (*skyCoverage_p)-=(*skyCoverage_p);
215 : // }
216 491 : sumWeight.resize(npol, nchan);
217 :
218 491 : convSupport=0;
219 :
220 491 : uvScale.resize(2);
221 491 : uvScale=0.0;
222 491 : uvScale(0)=Float(nx)*image->coordinates().increment()(0);
223 491 : uvScale(1)=Float(ny)*image->coordinates().increment()(1);
224 :
225 491 : uvOffset.resize(2);
226 491 : uvOffset(0)=nx/2;
227 491 : uvOffset(1)=ny/2;
228 :
229 982 : gridder.reset(new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
230 491 : uvScale, uvOffset,
231 491 : "SF"));
232 :
233 : // Set up image cache needed for gridding.
234 491 : if(imageCache) delete imageCache; imageCache=0;
235 : /*
236 : if(isTiled) {
237 : Float tileOverlap=0.5;
238 : tilesize=min(256,tilesize);
239 : IPosition tileShape=IPosition(4,tilesize,tilesize,npol,nchan);
240 : Vector<Float> tileOverlapVec(4);
241 : tileOverlapVec=0.0;
242 : tileOverlapVec(0)=tileOverlap;
243 : tileOverlapVec(1)=tileOverlap;
244 : Int tmpCacheVal=static_cast<Int>(cachesize);
245 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
246 : tileOverlapVec,
247 : (tileOverlap>0.0));
248 :
249 : }
250 : */
251 491 : }
252 :
253 : // This is nasty, we should use CountedPointers here.
254 628 : MosaicFT::~MosaicFT() {
255 628 : if(imageCache) delete imageCache; imageCache=0;
256 : // if(arrayLattice) delete arrayLattice; arrayLattice=0;
257 628 : }
258 :
259 :
260 380 : void MosaicFT::setConvFunc(CountedPtr<SimplePBConvFunc>& pbconvFunc){
261 :
262 :
263 380 : pbConvFunc_p=pbconvFunc;
264 :
265 :
266 380 : }
267 :
268 72 : CountedPtr<SimplePBConvFunc>& MosaicFT::getConvFunc(){
269 72 : return pbConvFunc_p;
270 : }
271 :
272 88349 : void MosaicFT::findConvFunction(const ImageInterface<Complex>& iimage,
273 : const vi::VisBuffer2& vb, const Matrix<Double>& /*rotateduvw*/, const bool /*ispsf*/) {
274 :
275 :
276 : //oversample if image is small
277 : //But not more than 5000 pixels
278 88349 : convSampling=(max(nx, ny) < 50) ? 100: Int(ceil(5000.0/max(nx, ny)));
279 88349 : if(convSampling <1)
280 0 : convSampling=1;
281 88349 : if(pbConvFunc_p.null())
282 0 : pbConvFunc_p=new SimplePBConvFunc();
283 88349 : if(sj_p)
284 88349 : pbConvFunc_p->setSkyJones(sj_p.get());
285 : ////TEST for HetArray only for now
286 88349 : if(pbConvFunc_p->name()=="HetArrayConvFunc"){
287 74901 : if(convSampling <10)
288 30452 : convSampling=10;
289 74901 : AipsrcValue<Int>::find (convSampling, "mosaic.oversampling", 10);
290 : }
291 88349 : if((pbConvFunc_p->getVBUtil()).null()){
292 0 : if(vbutil_p.null()){
293 0 : vbutil_p=new VisBufferUtil(vb);
294 : }
295 0 : pbConvFunc_p->setVBUtil(vbutil_p);
296 : }
297 : //cerr << "NELEMS " << interpVisFreq_p.nelements() << " lsr " << lsrFreq_p.nelements() << endl;
298 88349 : pbConvFunc_p->findConvFunction(iimage, vb, convSampling, interpVisFreq_p, convFunc, weightConvFunc_p, convSizePlanes_p, convSupportPlanes_p,
299 88349 : convPolMap_p, convChanMap_p, convRowMap_p, (useConjConvFunc_p && !toVis_p), MVDirection(-(movingDirShift_p.getAngle())), fixMovingSource_p);
300 :
301 : // cerr << "MAX of convFunc " << max(abs(convFunc)) << endl;
302 : //For now only use one size and support
303 88349 : if(convSizePlanes_p.nelements() ==0)
304 0 : convSize=0;
305 : else
306 88349 : convSize=max(convSizePlanes_p);
307 88349 : if(convSupportPlanes_p.nelements() ==0)
308 0 : convSupport=0;
309 : else
310 88349 : convSupport=max(convSupportPlanes_p);
311 :
312 88349 : }
313 :
314 116 : void MosaicFT::initializeToVis(ImageInterface<Complex>& iimage,
315 : const vi::VisBuffer2& vb)
316 : {
317 116 : image=&iimage;
318 116 : toVis_p=true;
319 116 : ok();
320 :
321 : // if(convSize==0) {
322 116 : init(vb);
323 :
324 : // }
325 :
326 : // Initialize the maps for polarization and channel. These maps
327 : // translate visibility indices into image indices
328 116 : initMaps(vb);
329 116 : pbConvFunc_p->setVBUtil(vbutil_p);
330 116 : pbConvFunc_p->setUsePointing(usePointingTable_p);
331 : //make sure we rotate the first field too
332 116 : lastFieldId_p=-1;
333 116 : phaseShifter_p=new UVWMachine(*uvwMachine_p);
334 : //This is needed here as we need to know the grid correction before FFTing
335 116 : findConvFunction(*image, vb, vb.uvw());
336 :
337 116 : prepGridForDegrid();
338 :
339 116 : }
340 :
341 :
342 116 : void MosaicFT::prepGridForDegrid(){
343 :
344 : //For now isTiled=false
345 116 : isTiled=false;
346 116 : nx = image->shape()(0);
347 116 : ny = image->shape()(1);
348 116 : npol = image->shape()(2);
349 116 : nchan = image->shape()(3);
350 :
351 116 : IPosition gridShape(4, nx, ny, npol, nchan);
352 116 : griddedData.resize(gridShape);
353 116 : griddedData=Complex(0.0);
354 :
355 116 : IPosition stride(4, 1);
356 232 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
357 232 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
358 232 : IPosition trc(blc+image->shape()-stride);
359 :
360 116 : IPosition start(4, 0);
361 116 : griddedData(blc, trc) = image->getSlice(start, image->shape());
362 :
363 116 : image->clearCache();
364 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
365 116 : arrayLattice = new ArrayLattice<Complex>(griddedData);
366 116 : lattice=arrayLattice;
367 : {///UnDo the grid correction
368 116 : Int inx = lattice->shape()(0);
369 : //Int iny = lattice->shape()(1);
370 116 : Vector<Complex> correction(inx);
371 116 : correction=Complex(1.0, 0.0);
372 :
373 : // Int npixCorr= max(nx,ny);
374 116 : Vector<Float> sincConvX(nx);
375 78828 : for (Int ix=0;ix<nx;ix++) {
376 78712 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
377 78712 : if(ix==nx/2) {
378 116 : sincConvX(ix)=1.0;
379 : }
380 : else {
381 78596 : sincConvX(ix)=sin(x)/x;
382 : }
383 : }
384 116 : Vector<Float> sincConvY(ny);
385 78828 : for (Int ix=0;ix<ny;ix++) {
386 78712 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
387 78712 : if(ix==ny/2) {
388 116 : sincConvY(ix)=1.0;
389 : }
390 : else {
391 78596 : sincConvY(ix)=sin(x)/x;
392 : }
393 : }
394 :
395 116 : IPosition cursorShape(4, inx, 1, 1, 1);
396 116 : IPosition axisPath(4, 0, 1, 2, 3);
397 116 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
398 116 : LatticeIterator<Complex> lix(*lattice, lsx);
399 181980 : for(lix.reset();!lix.atEnd();lix++) {
400 : //Int pol=lix.position()(2);
401 : //Int chan=lix.position()(3);
402 :
403 181864 : Int iy=lix.position()(1);
404 : //gridder->correctX1D(correction,iy);
405 188082824 : for (Int ix=0;ix<nx;ix++) {
406 187900960 : correction(ix)=(sincConvX(ix)*sincConvY(iy));
407 : }
408 181864 : lix.rwVectorCursor()*=correction;
409 :
410 : }
411 116 : }
412 :
413 :
414 116 : logIO() << LogIO::DEBUGGING << "Starting FFT of image" << LogIO::POST;
415 : // Now do the FFT2D in place
416 116 : ft_p.c2cFFT(*lattice);
417 : ///////////////////////
418 : /*{
419 : CoordinateSystem ftCoords(image->coordinates());
420 : Int directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
421 : DirectionCoordinate dc=ftCoords.directionCoordinate(directionIndex);
422 : Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
423 : Vector<Int> shape(2); shape(0)=griddedData.shape()(0) ;shape(1)=griddedData.shape()(1);
424 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
425 : ftCoords.replaceCoordinate(*ftdc, directionIndex);
426 : delete ftdc; ftdc=0;
427 : PagedImage<Float> thisScreen(griddedData.shape(), ftCoords, String("MODEL_GRID_VIS"));
428 : thisScreen.put(amplitude(griddedData));
429 : }*/
430 : ////////////////////////
431 116 : logIO() << LogIO::DEBUGGING << "Finished FFT" << LogIO::POST;
432 :
433 :
434 :
435 116 : }
436 :
437 :
438 116 : void MosaicFT::finalizeToVis()
439 : {
440 116 : logIO() << LogOrigin("MosaicFT", "finalizeToVis") << LogIO::NORMAL;
441 116 : logIO()<< LogIO::NORMAL2 << "Time degrid " << timedegrid_p << LogIO::POST;
442 116 : timedegrid_p=0.0;
443 :
444 116 : if(!arrayLattice.null()) arrayLattice=0;
445 116 : if(!lattice.null()) lattice=0;
446 116 : griddedData.resize();
447 :
448 : /*
449 : if(isTiled) {
450 :
451 : logIO() << LogOrigin("MosaicFT", "finalizeToVis") << LogIO::NORMAL;
452 :
453 : AlwaysAssert(imageCache, AipsError);
454 : AlwaysAssert(image, AipsError);
455 : ostringstream o;
456 : imageCache->flush();
457 : imageCache->showCacheStatistics(o);
458 : logIO() << o.str() << LogIO::POST;
459 : }
460 : */
461 116 : if(pointingToImage)
462 0 : delete pointingToImage;
463 116 : pointingToImage=0;
464 116 : }
465 :
466 :
467 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
468 : // grid.
469 :
470 :
471 :
472 :
473 :
474 375 : void MosaicFT::initializeToSky(ImageInterface<Complex>& iimage,
475 : Matrix<Float>& weight,
476 : const vi::VisBuffer2& vb)
477 : {
478 : // image always points to the image
479 375 : image=&iimage;
480 375 : toVis_p=False;
481 : // if(convSize==0) {
482 375 : init(vb);
483 :
484 : // }
485 :
486 : // Initialize the maps for polarization and channel. These maps
487 : // translate visibility indices into image indices
488 375 : initMaps(vb);
489 375 : pbConvFunc_p->setVBUtil(vbutil_p);
490 375 : pbConvFunc_p->setUsePointing(usePointingTable_p);
491 : //make sure we rotate the first field too
492 375 : lastFieldId_p=-1;
493 375 : phaseShifter_p=new UVWMachine(*uvwMachine_p);
494 : //findConvFunction(*image, vb);
495 : /*if((image->shape().product())>cachesize) {
496 : isTiled=true;
497 : }
498 : else {
499 : isTiled=false;
500 : }
501 : */
502 : //For now isTiled has to be false
503 375 : isTiled=false;
504 375 : nx = image->shape()(0);
505 375 : ny = image->shape()(1);
506 375 : npol = image->shape()(2);
507 375 : nchan = image->shape()(3);
508 :
509 375 : sumWeight=0.0;
510 375 : weight.resize(sumWeight.shape());
511 375 : weight=0.0;
512 :
513 375 : image->clearCache();
514 : // Initialize for in memory or to disk gridding. lattice will
515 : // point to the appropriate Lattice, either the ArrayLattice for
516 : // in memory gridding or to the image for to disk gridding.
517 : /*if(isTiled) {
518 : imageCache->flush();
519 : image->set(Complex(0.0));
520 : lattice=CountedPtr<Lattice<Complex> >(image, false);
521 : if( !doneWeightImage_p && (convWeightImage_p==0)){
522 :
523 : convWeightImage_p=new TempImage<Complex> (iimage.shape(),
524 : iimage.coordinates());
525 :
526 :
527 :
528 :
529 : convWeightImage_p->set(Complex(0.0));
530 : weightLattice=convWeightImage_p;
531 :
532 : }
533 : }
534 : else*/
535 : {
536 375 : IPosition gridShape(4, nx, ny, npol, nchan);
537 375 : if(!useDoubleGrid_p) {
538 0 : griddedData.resize(gridShape);
539 0 : griddedData=Complex(0.0);
540 : }
541 : else {
542 375 : griddedData2.resize(gridShape);
543 375 : griddedData2=DComplex(0.0);
544 : }
545 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
546 : //arrayLattice = new ArrayLattice<Complex>(griddedData);
547 : //lattice=arrayLattice;
548 :
549 375 : if( !doneWeightImage_p && (convWeightImage_p==0)){
550 :
551 :
552 :
553 530 : convWeightImage_p=new TempImage<Complex> (iimage.shape(),
554 265 : iimage.coordinates());
555 265 : griddedWeight.resize(gridShape);
556 : /*IPosition stride(4, 1);
557 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
558 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
559 : IPosition trc(blc+image->shape()-stride);
560 :
561 : griddedWeight(blc, trc).set(Complex(0.0));
562 : */
563 265 : if(useDoubleGrid_p){
564 265 : griddedWeight2.resize(gridShape);
565 265 : griddedWeight2=DComplex(0.0);
566 : }
567 : else{
568 0 : griddedWeight=Complex(0.0);
569 : }
570 : //if(weightLattice) delete weightLattice; weightLattice=0;
571 265 : weightLattice = new ArrayLattice<Complex>(griddedWeight);
572 :
573 : }
574 :
575 375 : }
576 :
577 : //cerr << "initializetosky lastfield " << lastFieldId_p << endl;
578 : // AlwaysAssert(lattice, AipsError);
579 :
580 375 : }
581 :
582 0 : void MosaicFT::reset(){
583 : //call the base class reset
584 0 : FTMachine::reset();
585 0 : doneWeightImage_p=false;
586 0 : convWeightImage_p=nullptr;
587 0 : pbConvFunc_p->reset();
588 0 : }
589 :
590 375 : void MosaicFT::finalizeToSky()
591 : {
592 375 : logIO() << LogOrigin("MosaicFT", "finalizeToSky") << LogIO::NORMAL;
593 375 : logIO() << LogIO::NORMAL2 << "time to massage data " << timemass_p << LogIO::POST;
594 375 : logIO() << LogIO::NORMAL2<< "time gridding " << timegrid_p << LogIO::POST;
595 375 : timemass_p=0.0;
596 375 : timegrid_p=0.0;
597 : // Now we flush the cache and report statistics
598 : // For memory based, we don't write anything out yet.
599 : /*if(isTiled) {
600 : logIO() << LogOrigin("MosaicFT", "finalizeToSky") << LogIO::NORMAL;
601 :
602 : AlwaysAssert(image, AipsError);
603 : AlwaysAssert(imageCache, AipsError);
604 : imageCache->flush();
605 : ostringstream o;
606 : imageCache->showCacheStatistics(o);
607 : logIO() << o.str() << LogIO::POST;
608 : }
609 : */
610 :
611 :
612 375 : if(!doneWeightImage_p){
613 265 : if(useDoubleGrid_p){
614 265 : convertArray(griddedWeight, griddedWeight2);
615 : //Don't need the double-prec grid anymore...
616 265 : griddedWeight2.resize();
617 : }
618 265 : ft_p.c2cFFT(*weightLattice, false);
619 : //Get the stokes right
620 265 : CoordinateSystem coords=convWeightImage_p->coordinates();
621 265 : Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
622 265 : Int npol=1;
623 265 : Vector<Int> whichStokes(npol);
624 294 : if(stokes_p=="I" || stokes_p=="RR" || stokes_p=="LL" ||stokes_p=="XX"
625 294 : || stokes_p=="YY" || stokes_p=="Q" || stokes_p=="U" || stokes_p=="V"){
626 248 : npol=1;
627 248 : whichStokes(0)=Stokes::type(stokes_p);
628 : // if single plane Q U or V are used...the weight should be the I weight
629 : //if(stokes_p=="Q" || stokes_p=="U" || stokes_p=="V")
630 : //whichStokes(0)=Stokes::type("I");
631 : }
632 17 : else if(stokes_p=="IV"){
633 5 : npol=2;
634 5 : whichStokes.resize(2);
635 5 : whichStokes(0)=Stokes::I;
636 5 : whichStokes(1)=Stokes::V;
637 : }
638 12 : else if(stokes_p=="QU"){
639 4 : npol=2;
640 4 : whichStokes.resize(2);
641 4 : whichStokes(0)=Stokes::Q;
642 4 : whichStokes(1)=Stokes::U;
643 : }
644 8 : else if(stokes_p=="RRLL"){
645 0 : npol=2;
646 0 : whichStokes.resize(2);
647 0 : whichStokes(0)=Stokes::RR;
648 0 : whichStokes(1)=Stokes::LL;
649 : }
650 8 : else if(stokes_p=="XXYY"){
651 0 : npol=2;
652 0 : whichStokes.resize(2);
653 0 : whichStokes(0)=Stokes::XX;
654 0 : whichStokes(1)=Stokes::YY;
655 : }
656 8 : else if(stokes_p=="IQU"){
657 0 : npol=3;
658 0 : whichStokes.resize(3);
659 0 : whichStokes(0)=Stokes::I;
660 0 : whichStokes(1)=Stokes::Q;
661 0 : whichStokes(2)=Stokes::U;
662 : }
663 8 : else if(stokes_p=="IQUV"){
664 8 : npol=4;
665 8 : whichStokes.resize(4);
666 8 : whichStokes(0)=Stokes::I;
667 8 : whichStokes(1)=Stokes::Q;
668 8 : whichStokes(2)=Stokes::U;
669 8 : whichStokes(3)=Stokes::V;
670 : }
671 :
672 265 : StokesCoordinate newStokesCoord(whichStokes);
673 265 : coords.replaceCoordinate(newStokesCoord, stokesIndex);
674 265 : IPosition imshp=convWeightImage_p->shape();
675 265 : imshp(2)=npol;
676 :
677 :
678 265 : skyCoverage_p=new TempImage<Float> (imshp, coords,1.0);
679 530 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
680 530 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
681 265 : IPosition stride(4, 1);
682 530 : IPosition trc(blc+image->shape()-stride);
683 :
684 : // Do the copy
685 265 : IPosition start(4, 0);
686 265 : convWeightImage_p->put(griddedWeight(blc, trc));
687 265 : StokesImageUtil::ToStokesPSF(*skyCoverage_p, *convWeightImage_p);
688 265 : if(npol>1){
689 : // only the I get it right Q and U or V may end up with zero depending
690 : // if RR or XX
691 17 : blc(0)=0; blc(1)=0; blc(3)=0;blc(2)=0;
692 17 : trc=skyCoverage_p->shape()-stride;
693 17 : trc(2)=0;
694 34 : SubImage<Float> isubim(*skyCoverage_p, Slicer(blc, trc, Slicer::endIsLast));
695 50 : for (Int k=1; k < npol; ++k){
696 33 : blc(2)=k; trc(2)=k;
697 66 : SubImage<Float> quvsubim(*skyCoverage_p, Slicer(blc, trc, Slicer::endIsLast), true);
698 33 : quvsubim.copyData(isubim);
699 33 : }
700 :
701 17 : }
702 : //Store this image in the pbconvfunc object as
703 : //it can be used for rescaling or shared by other ftmachines that use
704 : //this pbconvfunc
705 265 : pbConvFunc_p->setWeightImage(skyCoverage_p);
706 265 : if(convWeightImage_p) delete convWeightImage_p;
707 265 : convWeightImage_p=0;
708 265 : doneWeightImage_p=true;
709 :
710 : /*
711 : if(0){
712 : PagedImage<Float> thisScreen(skyCoverage_p->shape(),
713 : skyCoverage_p->coordinates(), "Screen");
714 : thisScreen.copyData(*skyCoverage_p);
715 : }
716 : */
717 :
718 265 : }
719 :
720 375 : if(!weightLattice.null()) weightLattice=0;
721 375 : griddedWeight.resize();
722 375 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
723 375 : }
724 :
725 0 : void MosaicFT::setWeightImage(CountedPtr<ImageInterface<Float> >& wgtimage){
726 0 : IPosition shp=wgtimage->shape();
727 0 : CoordinateSystem cs=wgtimage->coordinates();
728 0 : CountedPtr<TempImage<Float> > wgtim=new TempImage<Float>(shp, cs);
729 0 : wgtim->copyData(*(wgtimage));
730 0 : skyCoverage_p=wgtim;
731 0 : Record rec=skyCoverage_p->miscInfo();
732 : //For mosaicFTNew it has the nx*ny factor already in
733 0 : rec.define("isscaled", True);
734 0 : skyCoverage_p->setMiscInfo(rec);
735 : //cerr << "IN SET " << max(wgtimage->get()) << endl;
736 0 : pbConvFunc_p->setWeightImage(skyCoverage_p);
737 0 : doneWeightImage_p=true;
738 0 : }
739 :
740 0 : Array<Complex>* MosaicFT::getDataPointer(const IPosition& centerLoc2D,
741 : Bool readonly) {
742 : Array<Complex>* result;
743 : // Is tiled: get tiles and set up offsets
744 0 : centerLoc(0)=centerLoc2D(0);
745 0 : centerLoc(1)=centerLoc2D(1);
746 0 : result=&imageCache->tile(offsetLoc,centerLoc, readonly);
747 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
748 0 : return result;
749 : }
750 :
751 : #define NEED_UNDERSCORES
752 : #if defined(NEED_UNDERSCORES)
753 : #define sectgmoss2 sectgmoss2_
754 : #define gmoss2 gmoss2_
755 : #define sectgmosd2 sectgmosd2_
756 : #define gmosd2 gmosd2_
757 : #define sectdmos2 sectdmos2_
758 : #define dmos2 dmos2_
759 : #define gmoswgtd gmoswgtd_
760 : #define gmoswgts gmoswgts_
761 : #define locuvw locuvw_
762 : #endif
763 :
764 : extern "C" {
765 : void locuvw(const Double*, const Double*, const Double*, const Int*, const Double*, const Double*, const Int*,
766 : Int*, Int*, Complex*, const Int*, const Int*, const Double*);
767 : void gmoswgtd(const Int*/*nvispol*/, const Int*/*nvischan*/,
768 : const Int*/*flag*/, const Int*/*rflag*/, const Float*/*weight*/, const Int*/*nrow*/,
769 : const Int*/*nx*/, const Int*/*ny*/, const Int*/*npol*/, const Int*/*nchan*/,
770 : const Int*/*support*/, const Int*/*convsize*/, const Int*/*sampling*/,
771 : const Int*/*chanmap*/, const Int*/*polmap*/,
772 : DComplex* /*weightgrid*/, Double* /*sumwt*/, const Complex*/*convweight*/, const Int*/*convplanemap*/,
773 : const Int*/*convchanmap*/, const Int*/*convpolmap*/,
774 : const Int*/*nconvplane*/, const Int*/*nconvchan*/, const Int*/*nconvpol*/, const Int*/*rbeg*/,
775 : const Int*/*rend*/, const Int*/*loc*/, const Int*/*off*/, const Complex*/*phasor*/);
776 : void gmoswgts(const Int*/*nvispol*/, const Int*/*nvischan*/,
777 : const Int*/*flag*/, const Int*/*rflag*/, const Float*/*weight*/, const Int*/*nrow*/,
778 : const Int*/*nx*/, const Int*/*ny*/, const Int*/*npol*/, const Int*/*nchan*/,
779 : const Int*/*support*/, const Int*/*convsize*/, const Int*/*sampling*/,
780 : const Int*/*chanmap*/, const Int*/*polmap*/,
781 : Complex* /*weightgrid*/, Double*/*sumwt*/, const Complex*/*convweight*/, const Int*/*convplanemap*/,
782 : const Int*/*convchanmap*/, const Int*/*convpolmap*/,
783 : const Int*/*nconvplane*/, const Int*/*nconvchan*/, const Int*/*nconvpol*/, const Int*/*rbeg*/,
784 : const Int*/*rend*/, const Int*/*loc*/, const Int*/*off*/, const Complex*/*phasor*/);
785 : void sectgmosd2(const Complex* /*values*/,
786 : Int* /*nvispol*/, Int* /*nvischan*/,
787 : Int* /*dopsf*/, const Int* /*flag*/, const Int* /*rflag*/, const Float* /*weight*/,
788 : Int* /* nrow*/, DComplex* /*grid*/, Int* /*nx*/, Int* /*ny*/, Int * /*npol*/, Int * /*nchan */,
789 : Int*/*support*/, Int*/*convsize*/, Int*/*sampling*/, const Complex*/*convfunc*/,
790 : const Int*/*chanmap*/, const Int*/*polmap*/,
791 : Double*/*sumwgt*/, const Int*/*convplanemap*/,
792 : const Int*/*convchanmap*/, const Int*/*convpolmap*/,
793 : Int*/*nconvplane*/, Int*/*nconvchan*/, Int* /*nconvpol*/,
794 : const Int*/*x0*/,const Int*/*y0*/, const Int*/*nxsub*/, const Int*/*nysub*/, const Int*/*rbeg*/,
795 : const Int* /*rend*/, const Int*/*loc*/, const Int* /*off*/, const Complex*/*phasor*/);
796 :
797 : void sectgmoss2(const Complex* /*values*/,
798 : Int* /*nvispol*/, Int* /*nvischan*/,
799 : Int* /*dopsf*/, const Int* /*flag*/, const Int* /*rflag*/, const Float* /*weight*/,
800 : Int* /* nrow*/, Complex* /*grid*/, Int* /*nx*/, Int* /*ny*/, Int * /*npol*/, Int * /*nchan */,
801 : Int*/*support*/, Int*/*convsize*/, Int*/*sampling*/, const Complex*/*convfunc*/,
802 : const Int*/*chanmap*/, const Int*/*polmap*/,
803 : Double*/*sumwgt*/, const Int*/*convplanemap*/,
804 : const Int*/*convchanmap*/, const Int*/*convpolmap*/,
805 : Int*/*nconvplane*/, Int*/*nconvchan*/, Int* /*nconvpol*/,
806 : const Int*/*x0*/,const Int*/*y0*/, const Int*/*nxsub*/, const Int*/*nysub*/, const Int*/*rbeg*/,
807 : const Int* /*rend*/, const Int*/*loc*/, const Int* /*off*/, const Complex*/*phasor*/);
808 :
809 :
810 : void gmosd2(const Double*,
811 : Double*,
812 : const Complex*,
813 : Int*,
814 : Int*,
815 : Int*,
816 : const Int*,
817 : const Int*,
818 : const Float*,
819 : Int*,
820 : Int*,
821 : Double*,
822 : Double*,
823 : DComplex*,
824 : Int*,
825 : Int*,
826 : Int *,
827 : Int *,
828 : const Double*,
829 : const Double*,
830 : Int*,
831 : Int*,
832 : Int*,
833 : const Complex*,
834 : Int*,
835 : Int*,
836 : Double*,
837 : DComplex*,
838 : Complex*,
839 : Int*,
840 : Int*,
841 : Int*,Int*, Int*, Int*, Int*);
842 : /* void gmoss(const Double*,
843 : Double*,
844 : const Complex*,
845 : Int*,
846 : Int*,
847 : Int*,
848 : const Int*,
849 : const Int*,
850 : const Float*,
851 : Int*,
852 : Int*,
853 : Double*,
854 : Double*,
855 : Complex*,
856 : Int*,
857 : Int*,
858 : Int *,
859 : Int *,
860 : const Double*,
861 : const Double*,
862 : Int*,
863 : Int*,
864 : Int*,
865 : const Complex*,
866 : Int*,
867 : Int*,
868 : Double*,
869 : Complex*,
870 : Complex*,
871 : Int*,
872 : Int*,
873 : Int*);
874 : */
875 : void gmoss2(const Double*,
876 : Double*,
877 : const Complex*,
878 : Int*,
879 : Int*,
880 : Int*,
881 : const Int*,
882 : const Int*,
883 : const Float*,
884 : Int*, //10
885 : Int*,
886 : Double*,
887 : Double*,
888 : Complex*,
889 : Int*,
890 : Int*,
891 : Int *,
892 : Int *,
893 : const Double*,
894 : const Double*, //20
895 : Int*,
896 : Int*,
897 : Int*,
898 : const Complex*,
899 : Int*,
900 : Int*,
901 : Double*,
902 : Complex*,
903 : Complex*,
904 : Int*, //30
905 : Int*,
906 : Int*, Int*, Int*, Int*, Int*);
907 :
908 : /* void dmos(const Double*,
909 : Double*,
910 : Complex*,
911 : Int*,
912 : Int*,
913 : const Int*,
914 : const Int*,
915 : Int*,
916 : Int*,
917 : Double*, //10
918 : Double*,
919 : const Complex*,
920 : Int*,
921 : Int*,
922 : Int *,
923 : Int *,
924 : const Double*,
925 : const Double*,
926 : Int*,
927 : Int*,//20
928 : Int*,
929 : const Complex*,
930 : Int*,
931 : Int*,
932 : Int*,
933 : Int*);
934 : */
935 :
936 : void dmos2(const Double*,
937 : Double*,
938 : Complex*,
939 : Int*,
940 : Int*,
941 : const Int*,
942 : const Int*,
943 : Int*,
944 : Int*,
945 : Double*,
946 : Double*,
947 : const Complex*,
948 : Int*,
949 : Int*,
950 : Int *,
951 : Int *,
952 : const Double*,
953 : const Double*,
954 : Int*,
955 : Int*,
956 : Int*,
957 : const Complex*,
958 : Int*,
959 : Int*,
960 : Int*,
961 : Int*, Int*, Int*, Int*, Int*);
962 : void sectdmos2(Complex*,
963 : Int*,
964 : Int*,
965 : const Int*,
966 : const Int*,
967 : Int*,
968 : const Complex*,
969 : Int*,
970 : Int*,
971 : Int *,
972 : Int *,
973 : Int*,
974 : Int*,
975 : Int*,
976 : const Complex*,
977 : const Int*,
978 : const Int*,
979 : const Int*,
980 : const Int*,
981 : const Int*,
982 : Int*, Int*, Int*,
983 : //rbeg
984 : const Int*,
985 : const Int*,
986 : const Int*,
987 : const Int*,
988 : const Complex*);
989 :
990 :
991 :
992 : }
993 76393 : void MosaicFT::put(const vi::VisBuffer2& vb, Int row, Bool dopsf,
994 : FTMachine::Type type)
995 : {
996 :
997 :
998 :
999 :
1000 76393 : Timer tim;
1001 76393 : tim.mark();
1002 :
1003 76393 : matchChannel(vb);
1004 :
1005 :
1006 : //cerr << "CHANMAP " << chanMap << endl;
1007 : //No point in reading data if its not matching in frequency
1008 76393 : if(max(chanMap)==-1)
1009 0 : return;
1010 :
1011 : //const Matrix<Float> *imagingweight;
1012 : //imagingweight=&(vb.imagingWeight());
1013 76393 : Matrix<Float> imagingweight;
1014 76393 : getImagingWeight(imagingweight, vb);
1015 :
1016 76393 : if(dopsf) type=FTMachine::PSF;
1017 :
1018 76393 : Cube<Complex> data;
1019 : //Fortran gridder need the flag as ints
1020 76393 : Cube<Int> flags;
1021 76393 : Matrix<Float> elWeight;
1022 76393 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
1023 :
1024 :
1025 :
1026 : Bool iswgtCopy;
1027 : const Float *wgtStorage;
1028 76393 : wgtStorage=elWeight.getStorage(iswgtCopy);
1029 :
1030 :
1031 :
1032 :
1033 : Bool isCopy;
1034 76393 : const Complex *datStorage=0;
1035 :
1036 : // cerr << "dopsf " << dopsf << " isWeightCopy " << iswgtCopy << " " << wgtStorage<< endl;
1037 76393 : if(!dopsf)
1038 46760 : datStorage=data.getStorage(isCopy);
1039 :
1040 :
1041 : // If row is -1 then we pass through all rows
1042 : Int startRow, endRow, nRow;
1043 76393 : if (row==-1) {
1044 76393 : nRow=vb.nRows();
1045 76393 : startRow=0;
1046 76393 : endRow=nRow-1;
1047 : } else {
1048 0 : nRow=1;
1049 0 : startRow=row;
1050 0 : endRow=row;
1051 : }
1052 :
1053 : // Get the uvws in a form that Fortran can use and do that
1054 : // necessary phase rotation.
1055 76393 : Matrix<Double> uvw(negateUV(vb));
1056 76393 : Vector<Double> dphase(vb.nRows());
1057 76393 : dphase=0.0;
1058 :
1059 76393 : doUVWRotation_p=true;
1060 76393 : girarUVW(uvw, dphase, vb);
1061 76393 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1062 : // This needs to be after the interp to get the interpolated channels
1063 : //Also has to be after rotateuvw in case tracking is on
1064 : //cerr << "orig " << vb.uvw().row(2) << endl;
1065 : //vi::VisBuffer2& vbRotuvw=const_cast<vi::VisBuffer2&>(vb);
1066 : //vbRotuvw.setUvw(uvw);
1067 :
1068 : // cerr << "rot " << vbRotuvw.uvw().row(2) << endl;
1069 :
1070 76393 : findConvFunction(*image, vb, uvw, dopsf);
1071 :
1072 : //cerr << "Put convsup " << convSupport << " max min convFunc " << max(convFunc) << " " << min(convFunc) << " " << max(weightConvFunc_p) << min(weightConvFunc_p) << "SHP " << convFunc.shape() << " " << weightConvFunc_p.shape() << endl;
1073 : //cerr << "convRowMap " << convRowMap_p << " " << convChanMap_p << " " << convPolMap_p << endl;
1074 : //nothing to grid here as the pointing resulted in a zero support convfunc
1075 76393 : if(convSupport <= 0)
1076 0 : return;
1077 :
1078 : // Get the pointing positions. This can easily consume a lot
1079 : // of time thus we are for now assuming a field per
1080 : // vb chunk...need to change that accordingly if we start using
1081 : // multiple pointings per vb.
1082 : //Warning
1083 :
1084 : // Take care of translation of Bools to Integer
1085 76393 : Int idopsf=0;
1086 76393 : if(dopsf) idopsf=1;
1087 :
1088 :
1089 76393 : Vector<Int> rowFlags(vb.nRows());
1090 76393 : rowFlags=0;
1091 76393 : rowFlags(vb.flagRow())=true;
1092 76393 : if(!usezero_p) {
1093 16128505 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1094 16052112 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
1095 : }
1096 : }
1097 :
1098 :
1099 :
1100 : //cerr << "convSamp " << convSampling << " convsupp " << convSupport << " consize " << convSize << " convFunc " << convFunc.shape() << endl;
1101 : //TESTOO
1102 : /*{
1103 : ArrayIterator<Complex> itC(convFunc, IPosition(2,0,1));
1104 : ArrayIterator<Complex> itW(weightConvFunc_p, IPosition(2,0,1));
1105 : itC.origin();
1106 : itW.origin();
1107 : Int k=0;
1108 : while(!itC.pastEnd()){
1109 : cerr << k << "sum conv plane " << sum(itC.array()) << " wt " << sum(itW.array()) << endl;
1110 :
1111 : itC.next();
1112 : itW.next();
1113 : ++k;
1114 : }
1115 :
1116 : }*/
1117 : //TESTOO
1118 :
1119 : //Tell the gridder to grid the weights too ...need to do that once only
1120 : //Int doWeightGridding=1;
1121 : //if(doneWeightImage_p)
1122 : // doWeightGridding=-1;
1123 : Bool del;
1124 : // IPosition s(flags.shape());
1125 76393 : const IPosition& fs=flags.shape();
1126 : //cerr << "flags shape " << fs << endl;
1127 76393 : std::vector<Int>s(fs.begin(), fs.end());
1128 76393 : Int nvp=s[0];
1129 76393 : Int nvc=s[1];
1130 76393 : Int nvisrow=s[2];
1131 76393 : Int csamp=convSampling;
1132 : Bool uvwcopy;
1133 76393 : const Double *uvwstor=uvw.getStorage(uvwcopy);
1134 : Bool gridcopy;
1135 : Bool convcopy;
1136 : Bool wconvcopy;
1137 76393 : const Complex *convstor=convFunc.getStorage(convcopy);
1138 76393 : const Complex *wconvstor=weightConvFunc_p.getStorage(wconvcopy);
1139 76393 : Int nPolConv=convFunc.shape()[2];
1140 76393 : Int nChanConv=convFunc.shape()[3];
1141 76393 : Int nConvFunc=convFunc.shape()(4);
1142 : Bool weightcopy;
1143 : ////////**************************
1144 76393 : Cube<Int> loc(2, nvc, nRow);
1145 76393 : Cube<Int> off(2, nvc, nRow);
1146 76393 : Matrix<Complex> phasor(nvc, nRow);
1147 : Bool delphase;
1148 76393 : Complex * phasorstor=phasor.getStorage(delphase);
1149 76393 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1150 76393 : const Double * scalestor=uvScale.getStorage(del);
1151 76393 : const Double * offsetstor=uvOffset.getStorage(del);
1152 76393 : Int * locstor=loc.getStorage(del);
1153 76393 : Int * offstor=off.getStorage(del);
1154 76393 : const Double *dpstor=dphase.getStorage(del);
1155 : Int irow;
1156 76393 : Int nth=1;
1157 : #ifdef _OPENMP
1158 76393 : if(numthreads_p >0){
1159 0 : nth=min(numthreads_p, omp_get_max_threads());
1160 : }
1161 : else{
1162 76393 : nth= omp_get_max_threads();
1163 : }
1164 : //nth=min(4,nth);
1165 : #endif
1166 76393 : Double cinv=Double(1.0)/C::c;
1167 :
1168 76393 : Int dow=0;
1169 76393 : #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)
1170 : {
1171 : #pragma omp for
1172 : for (irow=startRow; irow<=endRow;irow++){
1173 : /*locateuvw(uvwstor,dpstor, visfreqstor, nvc, scalestor, offsetstor, csamp,
1174 : locstor,
1175 : offstor, phasorstor, irow, false);*/
1176 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
1177 : }
1178 :
1179 : }//end pragma parallel
1180 :
1181 :
1182 :
1183 76393 : timemass_p +=tim.real();
1184 : Int ixsub, iysub, icounter;
1185 76393 : ixsub=1;
1186 76393 : iysub=1;
1187 : //////***********************DEBUGGING
1188 : //nth=1;
1189 : ////////***************
1190 76393 : if (nth >3){
1191 76393 : ixsub=8;
1192 76393 : iysub=8;
1193 : }
1194 0 : else if(nth >1){
1195 0 : ixsub=2;
1196 0 : iysub=2;
1197 : }
1198 76393 : Int rbeg=startRow+1;
1199 76393 : Int rend=endRow+1;
1200 76393 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
1201 76393 : Vector<Double *> swgtptr(ixsub*iysub);
1202 76393 : Vector<Bool> swgtdel(ixsub*iysub);
1203 4965545 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1204 4889152 : sumwgt[icounter].resize(sumWeight.shape());
1205 4889152 : sumwgt[icounter].set(0.0);
1206 4889152 : swgtptr[icounter]=sumwgt[icounter].getStorage(swgtdel(icounter));
1207 : }
1208 : //cerr << "done thread " << doneThreadPartition_p << " " << ixsub*iysub << endl;
1209 76393 : if(doneThreadPartition_p < 0){
1210 265 : xsect_p.resize(ixsub*iysub);
1211 265 : ysect_p.resize(ixsub*iysub);
1212 265 : nxsect_p.resize(ixsub*iysub);
1213 265 : nysect_p.resize(ixsub*iysub);
1214 17225 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1215 16960 : findGridSector(nx, ny, ixsub, iysub, 0, 0, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
1216 : }
1217 : }
1218 76393 : Vector<Int> xsect, ysect, nxsect, nysect;
1219 76393 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
1220 : //cerr << xsect.shape() << " " << xsect << endl;
1221 76393 : const Int* pmapstor=polMap.getStorage(del);
1222 76393 : const Int* cmapstor=chanMap.getStorage(del);
1223 : // Dummy sumwt for gridweight part
1224 76393 : Matrix<Double> dumSumWeight(npol, nchan);
1225 76393 : dumSumWeight=sumWeight;
1226 : Bool isDSWC;
1227 76393 : Double *dsumwtstor=dumSumWeight.getStorage(isDSWC);
1228 76393 : Int nc=nchan;
1229 76393 : Int np=npol;
1230 76393 : Int nxp=nx;
1231 76393 : Int nyp=ny;
1232 76393 : Int csize=convSize;
1233 76393 : const Int * flagstor=flags.getStorage(del);
1234 76393 : const Int * rowflagstor=rowFlags.getStorage(del);
1235 76393 : Int csupp=convSupport;
1236 76393 : const Int *convrowmapstor=convRowMap_p.getStorage(del);
1237 76393 : const Int *convchanmapstor=convChanMap_p.getStorage(del);
1238 76393 : const Int *convpolmapstor=convPolMap_p.getStorage(del);
1239 : ///
1240 :
1241 :
1242 : ////////***************************
1243 76393 : tim.mark();
1244 :
1245 76393 : if(useDoubleGrid_p) {
1246 76393 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
1247 :
1248 76393 : #pragma omp parallel default(none) private(icounter, del) firstprivate(idopsf, /*doWeightGridding,*/ datStorage, wgtStorage, flagstor, rowflagstor, convstor, wconvstor, pmapstor, cmapstor, gridstor, csupp, nxp, nyp, np, nc,ixsub, iysub, rend, rbeg, csamp, csize, nvp, nvc, nvisrow, phasorstor, locstor, offstor, convrowmapstor, convchanmapstor, convpolmapstor, nPolConv, nChanConv, nConvFunc,xsect, ysect, nxsect, nysect) shared(swgtptr)
1249 : {
1250 : #pragma omp for schedule(dynamic)
1251 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
1252 : Int x0=xsect(icounter);
1253 : Int y0=ysect(icounter);
1254 : Int nxsub=nxsect(icounter);
1255 : Int nysub=nysect(icounter);
1256 :
1257 : /*
1258 : ix= (icounter+1)-((icounter)/ixsub)*ixsub;
1259 : iy=(icounter)/ixsub+1;
1260 : y0=(nyp/iysub)*(iy-1)+1;
1261 : nysub=nyp/iysub;
1262 : if( iy == iysub) {
1263 : nysub=nyp-(nyp/iysub)*(iy-1);
1264 : }
1265 : x0=(nxp/ixsub)*(ix-1)+1;
1266 : nxsub=nxp/ixsub;
1267 : if( ix == ixsub){
1268 : nxsub=nxp-(nxp/ixsub)*(ix-1);
1269 : }
1270 : */
1271 :
1272 : sectgmosd2(datStorage,
1273 : &nvp,
1274 : &nvc,
1275 : &idopsf,
1276 : flagstor,
1277 : rowflagstor,
1278 : wgtStorage,
1279 : &nvisrow,
1280 : gridstor,
1281 : &nxp,
1282 : &nyp,
1283 : &np,
1284 : &nc,
1285 : &csupp,
1286 : &csize,
1287 : &csamp,
1288 : convstor,
1289 : cmapstor,
1290 : pmapstor,
1291 : swgtptr[icounter],
1292 : convrowmapstor,
1293 : convchanmapstor,
1294 : convpolmapstor,
1295 : &nConvFunc, &nChanConv, &nPolConv,
1296 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
1297 : phasorstor
1298 : );
1299 : }
1300 : }//end pragma parallel
1301 4965545 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1302 4889152 : sumwgt[icounter].putStorage(swgtptr[icounter],swgtdel[icounter]);
1303 4889152 : sumWeight=sumWeight+sumwgt[icounter];
1304 : }
1305 :
1306 : //cerr << "SUMWEIG " << sumWeight << endl;
1307 76393 : griddedData2.putStorage(gridstor, gridcopy);
1308 76393 : if(dopsf && (nth >4))
1309 10536 : tweakGridSector(nx, ny, ixsub, iysub);
1310 76393 : timegrid_p+=tim.real();
1311 76393 : tim.mark();
1312 76393 : if(!doneWeightImage_p){
1313 : //This can be parallelized by making copy of the central part of the griddedWeight
1314 : //and adding it after dooing the gridding
1315 50203 : DComplex *gridwgtstor=griddedWeight2.getStorage(weightcopy);
1316 50203 : gmoswgtd(&nvp, &nvc,flagstor, rowflagstor, wgtStorage, &nvisrow,
1317 : &nxp, &nyp, &np, &nc, &csupp, &csize, &csamp,
1318 : cmapstor, pmapstor,
1319 : gridwgtstor, dsumwtstor, wconvstor, convrowmapstor,
1320 : convchanmapstor, convpolmapstor,
1321 : &nConvFunc, &nChanConv, &nPolConv, &rbeg,
1322 : &rend, locstor, offstor, phasorstor);
1323 50203 : griddedWeight2.putStorage(gridwgtstor, weightcopy);
1324 :
1325 : }
1326 76393 : timemass_p+=tim.real();
1327 : }
1328 : else {
1329 : //cerr << "maps " << convChanMap_p << " " << chanMap << endl;
1330 : //cerr << "nchan " << nchan << " nchanconv " << nChanConv << endl;
1331 0 : Complex *gridstor=griddedData.getStorage(gridcopy);
1332 0 : #pragma omp parallel default(none) private(icounter, del) firstprivate(idopsf, /*doWeightGridding,*/ datStorage, wgtStorage, flagstor, rowflagstor, convstor, wconvstor, pmapstor, cmapstor, gridstor, csupp, nxp, nyp, np, nc,ixsub, iysub, rend, rbeg, csamp, csize, nvp, nvc, nvisrow, phasorstor, locstor, offstor, convrowmapstor, convchanmapstor, convpolmapstor, nPolConv, nChanConv, nConvFunc, xsect, ysect, nxsect, nysect) shared(swgtptr)
1333 : {
1334 : #pragma omp for schedule(dynamic)
1335 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
1336 : /*ix= (icounter+1)-((icounter)/ixsub)*ixsub;
1337 : iy=(icounter)/ixsub+1;
1338 : y0=(nyp/iysub)*(iy-1)+1;
1339 : nysub=nyp/iysub;
1340 : if( iy == iysub) {
1341 : nysub=nyp-(nyp/iysub)*(iy-1);
1342 : }
1343 : x0=(nxp/ixsub)*(ix-1)+1;
1344 : nxsub=nxp/ixsub;
1345 : if( ix == ixsub){
1346 : nxsub=nxp-(nxp/ixsub)*(ix-1);
1347 : }
1348 : */
1349 : Int x0=xsect(icounter);
1350 : Int y0=ysect(icounter);
1351 : Int nxsub=nxsect(icounter);
1352 : Int nysub=nysect(icounter);
1353 :
1354 : sectgmoss2(datStorage,
1355 : &nvp,
1356 : &nvc,
1357 : &idopsf,
1358 : flagstor,
1359 : rowflagstor,
1360 : wgtStorage,
1361 : &nvisrow,
1362 : gridstor,
1363 : &nxp,
1364 : &nyp,
1365 : &np,
1366 : &nc,
1367 : &csupp,
1368 : &csize,
1369 : &csamp,
1370 : convstor,
1371 : cmapstor,
1372 : pmapstor,
1373 : swgtptr[icounter],
1374 : convrowmapstor,
1375 : convchanmapstor,
1376 : convpolmapstor,
1377 : &nConvFunc, &nChanConv, &nPolConv,
1378 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
1379 : phasorstor
1380 : );
1381 :
1382 :
1383 : }
1384 : } //end pragma
1385 0 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1386 0 : sumwgt[icounter].putStorage(swgtptr[icounter],swgtdel[icounter]);
1387 0 : sumWeight=sumWeight+sumwgt[icounter];
1388 : }
1389 0 : griddedData.putStorage(gridstor, gridcopy);
1390 0 : if(dopsf && (nth > 4))
1391 0 : tweakGridSector(nx, ny, ixsub, iysub);
1392 0 : timegrid_p+=tim.real();
1393 0 : tim.mark();
1394 0 : if(!doneWeightImage_p){
1395 0 : Complex *gridwgtstor=griddedWeight.getStorage(weightcopy);
1396 0 : gmoswgts(&nvp, &nvc,flagstor, rowflagstor, wgtStorage, &nvisrow,
1397 : &nxp, &nyp, &np, &nc, &csupp, &csize, &csamp,
1398 : cmapstor, pmapstor,
1399 : gridwgtstor, dsumwtstor, wconvstor, convrowmapstor,
1400 : convchanmapstor, convpolmapstor,
1401 : &nConvFunc, &nChanConv, &nPolConv, &rbeg,
1402 : &rend, locstor, offstor, phasorstor);
1403 0 : griddedWeight.putStorage(gridwgtstor, weightcopy);
1404 :
1405 : }
1406 0 : timemass_p+=tim.real();
1407 : }
1408 76393 : convFunc.freeStorage(convstor, convcopy);
1409 76393 : weightConvFunc_p.freeStorage(wconvstor, wconvcopy);
1410 76393 : dumSumWeight.putStorage(dsumwtstor, isDSWC);
1411 : //cerr << "dumSumwe " << dumSumWeight << endl;
1412 76393 : uvw.freeStorage(uvwstor, uvwcopy);
1413 76393 : if(!dopsf)
1414 46760 : data.freeStorage(datStorage, isCopy);
1415 :
1416 76393 : elWeight.freeStorage(wgtStorage,iswgtCopy);
1417 :
1418 :
1419 :
1420 :
1421 76393 : }
1422 :
1423 0 : void MosaicFT::gridImgWeights(const vi::VisBuffer2& vb){
1424 :
1425 0 : if(doneWeightImage_p)
1426 0 : return;
1427 0 : matchChannel(vb);
1428 :
1429 :
1430 : //cerr << "CHANMAP " << chanMap << endl;
1431 : //No point in reading data if its not matching in frequency
1432 0 : if(max(chanMap)==-1)
1433 0 : return;
1434 :
1435 : Int startRow, endRow, nRow;
1436 0 : nRow=vb.nRows();
1437 0 : startRow=0;
1438 0 : endRow=nRow-1;
1439 :
1440 :
1441 : //const Matrix<Float> *imagingweight;
1442 : //imagingweight=&(vb.imagingWeight());
1443 0 : Matrix<Float> imagingweight;
1444 0 : getImagingWeight(imagingweight, vb);
1445 :
1446 :
1447 0 : Cube<Complex> data;
1448 : //Fortran gridder need the flag as ints
1449 0 : Cube<Int> flags;
1450 0 : Matrix<Float> elWeight;
1451 0 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, FTMachine::PSF);
1452 :
1453 :
1454 :
1455 : Bool iswgtCopy;
1456 : const Float *wgtStorage;
1457 0 : wgtStorage=elWeight.getStorage(iswgtCopy);
1458 : Bool issumWgtCopy;
1459 0 : Double* sumwgtstor=sumWeight.getStorage(issumWgtCopy);
1460 :
1461 :
1462 :
1463 : // Get the uvws in a form that Fortran can use and do that
1464 : // necessary phase rotation.
1465 0 : Matrix<Double> uvw(negateUV(vb));
1466 0 : Vector<Double> dphase(vb.nRows());
1467 0 : dphase=0.0;
1468 :
1469 0 : doUVWRotation_p=true;
1470 0 : girarUVW(uvw, dphase, vb);
1471 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1472 : // This needs to be after the interp to get the interpolated channels
1473 : //Also has to be after rotateuvw in case tracking is on
1474 : //vi::VisBuffer2& vbRotuvw=const_cast<vi::VisBuffer2&>(vb);
1475 : //vbRotuvw.setUvw(uvw);
1476 :
1477 0 : findConvFunction(*image, vb, uvw);
1478 : //nothing to grid here as the pointing resulted in a zero support convfunc
1479 0 : if(convSupport <= 0)
1480 0 : return;
1481 :
1482 : Bool del;
1483 :
1484 0 : const Int* pmapstor=polMap.getStorage(del);
1485 0 : const Int* cmapstor=chanMap.getStorage(del);
1486 :
1487 0 : Vector<Int> rowFlags(vb.nRows());
1488 0 : rowFlags=0;
1489 0 : rowFlags(vb.flagRow())=true;
1490 0 : if(!usezero_p) {
1491 0 : for (uInt rownr=0; rownr< vb.nRows(); rownr++) {
1492 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
1493 : }
1494 : }
1495 :
1496 : //Fortran indexing
1497 :
1498 0 : Int rbeg=1;
1499 0 : Int rend=vb.nRows();
1500 :
1501 0 : const Int * flagstor=flags.getStorage(del);
1502 0 : const Int * rowflagstor=rowFlags.getStorage(del);
1503 :
1504 0 : const Int *convrowmapstor=convRowMap_p.getStorage(del);
1505 0 : const Int *convchanmapstor=convChanMap_p.getStorage(del);
1506 0 : const Int *convpolmapstor=convPolMap_p.getStorage(del);
1507 :
1508 : //Tell the gridder to grid the weights too ...need to do that once only
1509 : //Int doWeightGridding=1;
1510 : //if(doneWeightImage_p)
1511 : // doWeightGridding=-1;
1512 : // IPosition s(flags.shape());
1513 0 : const IPosition& fs=flags.shape();
1514 : //cerr << "flags shape " << fs << endl;
1515 0 : std::vector<Int>s(fs.begin(), fs.end());
1516 0 : Int nvp=s[0];
1517 0 : Int nvc=s[1];
1518 0 : Int nvisrow=s[2];
1519 0 : Int csamp=convSampling;
1520 : Bool uvwcopy;
1521 0 : const Double *uvwstor=uvw.getStorage(uvwcopy);
1522 : Bool gridcopy;
1523 : Bool convcopy;
1524 : Bool wconvcopy;
1525 0 : const Complex *wconvstor=weightConvFunc_p.getStorage(wconvcopy);
1526 0 : Int nPolConv=convFunc.shape()[2];
1527 0 : Int nChanConv=convFunc.shape()[3];
1528 0 : Int nConvFunc=convFunc.shape()(4);
1529 : Bool weightcopy;
1530 : ////////**************************
1531 0 : Cube<Int> loc(2, nvc, vb.nRows());
1532 0 : Cube<Int> off(2, nvc, vb.nRows());
1533 0 : Matrix<Complex> phasor(nvc, vb.nRows());
1534 : Bool delphase;
1535 0 : Complex * phasorstor=phasor.getStorage(delphase);
1536 0 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1537 0 : const Double * scalestor=uvScale.getStorage(del);
1538 0 : const Double * offsetstor=uvOffset.getStorage(del);
1539 0 : Int * locstor=loc.getStorage(del);
1540 0 : Int * offstor=off.getStorage(del);
1541 0 : const Double *dpstor=dphase.getStorage(del);
1542 :
1543 : Int irow;
1544 0 : Int nth=1;
1545 : #ifdef _OPENMP
1546 0 : if(numthreads_p >0){
1547 0 : nth=min(numthreads_p, omp_get_max_threads());
1548 : }
1549 : else{
1550 0 : nth= omp_get_max_threads();
1551 : }
1552 : //nth=min(4,nth);
1553 : #endif
1554 :
1555 0 : Double cinv=Double(1.0)/C::c;
1556 :
1557 0 : Int dow=0;
1558 :
1559 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)
1560 : {
1561 : #pragma omp for
1562 : for (irow=startRow; irow<=endRow;irow++){
1563 : /*locateuvw(uvwstor,dpstor, visfreqstor, nvc, scalestor, offsetstor, csamp,
1564 : locstor,
1565 : offstor, phasorstor, irow, false);*/
1566 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
1567 : }
1568 :
1569 : }//end pragma parallel
1570 :
1571 :
1572 :
1573 :
1574 0 : if(useDoubleGrid_p) {
1575 : //This can be parallelized by making copy of the central part of the griddedWeight
1576 : //and adding it after dooing the gridding
1577 0 : DComplex *gridwgtstor=griddedWeight2.getStorage(weightcopy);
1578 0 : gmoswgtd(&nvp, &nvc,flagstor, rowflagstor, wgtStorage, &nvisrow,
1579 0 : &nx, &ny, &npol, &nchan, &convSupport, &convSize, &convSampling,
1580 : cmapstor, pmapstor,
1581 : gridwgtstor, sumwgtstor, wconvstor, convrowmapstor,
1582 : convchanmapstor, convpolmapstor,
1583 : &nConvFunc, &nChanConv, &nPolConv, &rbeg,
1584 : &rend, locstor, offstor, phasorstor);
1585 0 : griddedWeight2.putStorage(gridwgtstor, weightcopy);
1586 :
1587 :
1588 :
1589 :
1590 :
1591 : }
1592 : else{
1593 0 : Complex *gridwgtstor=griddedWeight.getStorage(weightcopy);
1594 0 : gmoswgts(&nvp, &nvc,flagstor, rowflagstor, wgtStorage, &nvisrow,
1595 0 : &nx, &ny, &npol, &nchan, &convSupport, &convSize, &convSampling,
1596 : cmapstor, pmapstor,
1597 : gridwgtstor, sumwgtstor, wconvstor, convrowmapstor,
1598 : convchanmapstor, convpolmapstor,
1599 : &nConvFunc, &nChanConv, &nPolConv, &rbeg,
1600 : &rend, locstor, offstor, phasorstor);
1601 0 : griddedWeight.putStorage(gridwgtstor, weightcopy);
1602 :
1603 :
1604 : }
1605 0 : sumWeight.putStorage(sumwgtstor, issumWgtCopy);
1606 0 : elWeight.freeStorage(wgtStorage,iswgtCopy);
1607 :
1608 0 : }
1609 :
1610 22977 : void MosaicFT::get(vi::VisBuffer2& vb, Int row)
1611 : {
1612 :
1613 :
1614 :
1615 : // If row is -1 then we pass through all rows
1616 : Int startRow, endRow, nRow;
1617 22977 : if (row==-1) {
1618 22977 : nRow=vb.nRows();
1619 22977 : startRow=0;
1620 22977 : endRow=nRow-1;
1621 : // vb.modelVisCube()=Complex(0.0,0.0);
1622 : } else {
1623 0 : nRow=1;
1624 0 : startRow=row;
1625 0 : endRow=row;
1626 : // vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
1627 : }
1628 :
1629 :
1630 :
1631 :
1632 22977 : matchChannel(vb);
1633 :
1634 : //No point in reading data if its not matching in frequency
1635 22977 : if(max(chanMap)==-1)
1636 0 : return;
1637 :
1638 : // Get the uvws in a form that Fortran can use
1639 22977 : Matrix<Double> uvw(negateUV(vb));
1640 22977 : Vector<Double> dphase(vb.nRows());
1641 22977 : dphase=0.0;
1642 :
1643 22977 : doUVWRotation_p=true;
1644 22977 : girarUVW(uvw, dphase, vb);
1645 22977 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1646 :
1647 :
1648 :
1649 :
1650 22977 : Cube<Complex> data;
1651 22977 : Cube<Int> flags;
1652 22977 : getInterpolateArrays(vb, data, flags);
1653 :
1654 : //vi::VisBuffer2& vbRotuvw=const_cast<vi::VisBuffer2&>(vb);
1655 : //vbRotuvw.setUvw(uvw);
1656 :
1657 : //Need to get interpolated freqs
1658 22977 : findConvFunction(*image, vb, uvw);
1659 :
1660 : // no valid pointing in this buffer
1661 22977 : if(convSupport <= 0)
1662 0 : return;
1663 : Complex *datStorage;
1664 : Bool isCopy;
1665 22977 : datStorage=data.getStorage(isCopy);
1666 :
1667 :
1668 22977 : Vector<Int> rowFlags(vb.nRows());
1669 22977 : rowFlags=0;
1670 22977 : rowFlags(vb.flagRow())=true;
1671 22977 : if(!usezero_p) {
1672 3116103 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1673 3093126 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
1674 : }
1675 : }
1676 22977 : Int nvp=data.shape()[0];
1677 22977 : Int nvc=data.shape()[1];
1678 22977 : Int nvisrow=data.shape()[2];
1679 22977 : Int csamp=convSampling;
1680 22977 : Int csize=convSize;
1681 22977 : Int csupp=convSupport;
1682 22977 : Int nc=nchan;
1683 22977 : Int np=npol;
1684 22977 : Int nxp=nx;
1685 22977 : Int nyp=ny;
1686 : Bool uvwcopy;
1687 22977 : const Double *uvwstor=uvw.getStorage(uvwcopy);
1688 22977 : Int nPolConv=convFunc.shape()[2];
1689 22977 : Int nChanConv=convFunc.shape()[3];
1690 22977 : Int nConvFunc=convFunc.shape()(4);
1691 : ////////**************************
1692 22977 : Cube<Int> loc(2, nvc, nRow);
1693 22977 : Cube<Int> off(2, nvc, nRow);
1694 22977 : Matrix<Complex> phasor(nvc, nRow);
1695 : Bool delphase;
1696 : Bool del;
1697 22977 : const Int* pmapstor=polMap.getStorage(del);
1698 22977 : const Int* cmapstor=chanMap.getStorage(del);
1699 22977 : Complex * phasorstor=phasor.getStorage(delphase);
1700 22977 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1701 22977 : const Double * scalestor=uvScale.getStorage(del);
1702 22977 : const Double * offsetstor=uvOffset.getStorage(del);
1703 22977 : const Int * flagstor=flags.getStorage(del);
1704 22977 : const Int * rowflagstor=rowFlags.getStorage(del);
1705 22977 : Int * locstor=loc.getStorage(del);
1706 22977 : Int * offstor=off.getStorage(del);
1707 22977 : const Double *dpstor=dphase.getStorage(del);
1708 22977 : const Int *convrowmapstor=convRowMap_p.getStorage(del);
1709 22977 : const Int *convchanmapstor=convChanMap_p.getStorage(del);
1710 22977 : const Int *convpolmapstor=convPolMap_p.getStorage(del);
1711 : ////////***************************
1712 :
1713 : Int irow;
1714 22977 : Int nth=1;
1715 : #ifdef _OPENMP
1716 22977 : if(numthreads_p >0){
1717 0 : nth=min(numthreads_p, omp_get_max_threads());
1718 : }
1719 : else{
1720 22977 : nth= omp_get_max_threads();
1721 : }
1722 : //nth=min(4,nth);
1723 : #endif
1724 :
1725 22977 : Timer tim;
1726 22977 : tim.mark();
1727 :
1728 22977 : Int dow=0;
1729 22977 : Double cinv=Double(1.0)/C::c;
1730 22977 : #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)
1731 : {
1732 : #pragma omp for
1733 : for (irow=startRow; irow<=endRow;irow++){
1734 : /////////////////*locateuvw(uvwstor,dpstor, visfreqstor, nvc, scalestor, offsetstor, csamp,
1735 : // locstor,
1736 : /////////// offstor, phasorstor, irow, false);
1737 : //using the fortran version which is significantly faster ...this can account for 10% less overall degridding time
1738 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor,
1739 : offstor, phasorstor, &irow, &dow, &cinv);
1740 : }
1741 :
1742 : }//end pragma parallel
1743 22977 : Int rbeg=startRow+1;
1744 22977 : Int rend=endRow+1;
1745 22977 : Int npart=nth;
1746 :
1747 : Bool gridcopy;
1748 22977 : const Complex *gridstor=griddedData.getStorage(gridcopy);
1749 : Bool convcopy;
1750 : ////Degridding needs the conjugate ...doing it here
1751 22977 : Array<Complex> conjConvFunc=conj(convFunc);
1752 22977 : const Complex *convstor=conjConvFunc.getStorage(convcopy);
1753 22977 : Int ix=0;
1754 22977 : #pragma omp parallel default(none) private(ix, rbeg, rend) firstprivate(uvwstor, datStorage, flagstor, rowflagstor, convstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc, csamp, csize, csupp, nvp, nvc, nvisrow, phasorstor, locstor, offstor, nPolConv, nChanConv, nConvFunc, convrowmapstor, convpolmapstor, convchanmapstor, npart) num_threads(npart)
1755 : {
1756 : #pragma omp for schedule(dynamic)
1757 : for (ix=0; ix< npart; ++ix){
1758 : rbeg=ix*(nvisrow/npart)+1;
1759 : rend=(ix != (npart-1)) ? (rbeg+(nvisrow/npart)-1) : (rbeg+(nvisrow/npart)-1+nvisrow%npart) ;
1760 : //cerr << "maps " << convChanMap_p << " " << chanMap << endl;
1761 : //cerr << "nchan " << nchan << " nchanconv " << nChanConv << " npolconv " << nPolConv << " nRowConv " << nConvFunc << endl;
1762 : sectdmos2(
1763 : datStorage,
1764 : &nvp,
1765 : &nvc,
1766 : flagstor,
1767 : rowflagstor,
1768 : &nvisrow,
1769 : gridstor,
1770 : &nxp,
1771 : &nyp,
1772 : &np,
1773 : &nc,
1774 : &csupp,
1775 : &csize,
1776 : &csamp,
1777 : convstor,
1778 : cmapstor,
1779 : pmapstor,
1780 : convrowmapstor, convchanmapstor,
1781 : convpolmapstor,
1782 : &nConvFunc, &nChanConv, &nPolConv,
1783 : &rbeg, &rend, locstor, offstor, phasorstor
1784 : );
1785 :
1786 :
1787 : }
1788 : }//end pragma omp
1789 :
1790 :
1791 22977 : data.putStorage(datStorage, isCopy);
1792 22977 : griddedData.freeStorage(gridstor, gridcopy);
1793 22977 : convFunc.freeStorage(convstor, convcopy);
1794 :
1795 22977 : timedegrid_p+=tim.real();
1796 :
1797 22977 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1798 22977 : }
1799 :
1800 :
1801 : /*
1802 : void MosaicFT::get(VisBuffer& vb, Int row)
1803 : {
1804 :
1805 :
1806 :
1807 : // If row is -1 then we pass through all rows
1808 : Int startRow, endRow, nRow;
1809 : if (row==-1) {
1810 : nRow=vb.nRow();
1811 : startRow=0;
1812 : endRow=nRow-1;
1813 : // vb.modelVisCube()=Complex(0.0,0.0);
1814 : } else {
1815 : nRow=1;
1816 : startRow=row;
1817 : endRow=row;
1818 : // vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
1819 : }
1820 :
1821 :
1822 :
1823 :
1824 :
1825 : // Get the uvws in a form that Fortran can use
1826 : Matrix<Double> uvw(3, vb.uvw().nelements());
1827 : uvw=0.0;
1828 : Vector<Double> dphase(vb.uvw().nelements());
1829 : dphase=0.0;
1830 : //NEGATING to correct for an image inversion problem
1831 : for (Int i=startRow;i<=endRow;i++) {
1832 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
1833 : uvw(2,i)=vb.uvw()(i)(2);
1834 : }
1835 :
1836 : doUVWRotation_p=true;
1837 : girarUVW(uvw, dphase, vb);
1838 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1839 :
1840 :
1841 : //Check if ms has changed then cache new spw and chan selection
1842 : if(vb.newMS())
1843 : matchAllSpwChans(vb);
1844 :
1845 : //Here we redo the match or use previous match
1846 :
1847 : //Channel matching for the actual spectral window of buffer
1848 : if(doConversion_p[vb.spectralWindow()]){
1849 : matchChannel(vb.spectralWindow(), vb);
1850 : }
1851 : else{
1852 : chanMap.resize();
1853 : chanMap=multiChanMap_p[vb.spectralWindow()];
1854 : }
1855 : //No point in reading data if its not matching in frequency
1856 : if(max(chanMap)==-1)
1857 : return;
1858 :
1859 : Cube<Complex> data;
1860 : Cube<Int> flags;
1861 : getInterpolateArrays(vb, data, flags);
1862 : //Need to get interpolated freqs
1863 : findConvFunction(*image, vb);
1864 : // no valid pointing in this buffer
1865 : if(convSupport <= 0)
1866 : return;
1867 : Complex *datStorage;
1868 : Bool isCopy;
1869 : datStorage=data.getStorage(isCopy);
1870 :
1871 :
1872 : Vector<Int> rowFlags(vb.nRow());
1873 : rowFlags=0;
1874 : rowFlags(vb.flagRow())=true;
1875 : if(!usezero_p) {
1876 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1877 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
1878 : }
1879 : }
1880 : Int nPolConv=convFunc.shape()[2];
1881 : Int nChanConv=convFunc.shape()[3];
1882 : Int nConvFunc=convFunc.shape()(4);
1883 :
1884 : {
1885 : Bool del;
1886 : Bool uvwcopy;
1887 : const Double *uvwstor=uvw.getStorage(uvwcopy);
1888 : Bool gridcopy;
1889 : const Complex *gridstor=griddedData.getStorage(gridcopy);
1890 : Bool convcopy;
1891 : ////Degridding needs the conjugate ...doing it here
1892 : Array<Complex> conjConvFunc=conj(convFunc);
1893 : const Complex *convstor=conjConvFunc.getStorage(convcopy);
1894 : // IPosition s(data.shape());
1895 : const IPosition& fs=data.shape();
1896 : std::vector<Int> s(fs.begin(), fs.end());
1897 : //cerr << "maps " << convChanMap_p << " " << chanMap << endl;
1898 : //cerr << "nchan " << nchan << " nchanconv " << nChanConv << " npolconv " << nPolConv << " nRowConv " << nConvFunc << endl;
1899 : dmos2(uvwstor,
1900 : dphase.getStorage(del),
1901 : datStorage,
1902 : &s[0],
1903 : &s[1],
1904 : flags.getStorage(del),
1905 : rowFlags.getStorage(del),
1906 : &s[2],
1907 : &row,
1908 : uvScale.getStorage(del), //10
1909 : uvOffset.getStorage(del),
1910 : gridstor,
1911 : &nx,
1912 : &ny,
1913 : &npol,
1914 : &nchan,
1915 : interpVisFreq_p.getStorage(del),
1916 : &C::c,
1917 : &convSupport,
1918 : &convSize, //20
1919 : &convSampling,
1920 : convstor,
1921 : chanMap.getStorage(del),
1922 : polMap.getStorage(del),
1923 : convRowMap_p.getStorage(del), convChanMap_p.getStorage(del),
1924 : convPolMap_p.getStorage(del),
1925 : &nConvFunc, &nChanConv, &nPolConv);
1926 : data.putStorage(datStorage, isCopy);
1927 : uvw.freeStorage(uvwstor, uvwcopy);
1928 : griddedData.freeStorage(gridstor, gridcopy);
1929 : convFunc.freeStorage(convstor, convcopy);
1930 : }
1931 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1932 : }
1933 :
1934 : */
1935 :
1936 :
1937 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1938 : // return the resulting image
1939 0 : ImageInterface<Complex>& MosaicFT::getImage(Matrix<Float>& weights,
1940 : Bool normalize)
1941 : {
1942 : //AlwaysAssert(lattice, AipsError);
1943 0 : AlwaysAssert(image, AipsError);
1944 :
1945 0 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
1946 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
1947 : }
1948 :
1949 0 : logIO() << LogOrigin("MosaicFT", "getImage") << LogIO::NORMAL;
1950 :
1951 0 : weights.resize(sumWeight.shape());
1952 0 : convertArray(weights, sumWeight);
1953 0 : SynthesisUtilMethods::getResource("mem peak in getImage");
1954 :
1955 : // If the weights are all zero then we cannot normalize
1956 : // otherwise we don't care.
1957 0 : if(max(weights)==0.0) {
1958 0 : if(normalize) {
1959 : logIO() << LogIO::SEVERE << "No useful data in MosaicFT: weights all zero"
1960 0 : << LogIO::POST;
1961 : }
1962 : else {
1963 : logIO() << LogIO::WARN << "No useful data in MosaicFT: weights all zero"
1964 0 : << LogIO::POST;
1965 : }
1966 : }
1967 : else {
1968 :
1969 : //const IPosition latticeShape = lattice->shape();
1970 :
1971 : logIO() << LogIO::DEBUGGING
1972 0 : << "Starting FFT and scaling of image" << LogIO::POST;
1973 0 : if(useDoubleGrid_p){
1974 0 : ArrayLattice<DComplex> darrayLattice(griddedData2);
1975 0 : ft_p.c2cFFT(darrayLattice,false);
1976 0 : griddedData.resize(griddedData2.shape());
1977 0 : convertArray(griddedData, griddedData2);
1978 :
1979 : //Don't need the double-prec grid anymore...
1980 0 : griddedData2.resize();
1981 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1982 0 : lattice=arrayLattice;
1983 :
1984 0 : }
1985 : else{
1986 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1987 0 : lattice=arrayLattice;
1988 0 : ft_p.c2cFFT(*lattice,false);
1989 : }
1990 : {////Do the grid correction
1991 0 : Int inx = lattice->shape()(0);
1992 : //Int iny = lattice->shape()(1);
1993 0 : Vector<Complex> correction(inx);
1994 0 : correction=Complex(1.0, 0.0);
1995 :
1996 : // Int npixCorr= max(nx,ny);
1997 0 : Vector<Float> sincConvX(nx);
1998 0 : for (Int ix=0;ix<nx;ix++) {
1999 0 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
2000 0 : if(ix==nx/2) {
2001 0 : sincConvX(ix)=1.0;
2002 : }
2003 : else {
2004 0 : sincConvX(ix)=sin(x)/x;
2005 : }
2006 : }
2007 0 : Vector<Float> sincConvY(ny);
2008 0 : for (Int ix=0;ix<ny;ix++) {
2009 0 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
2010 0 : if(ix==ny/2) {
2011 0 : sincConvY(ix)=1.0;
2012 : }
2013 : else {
2014 0 : sincConvY(ix)=sin(x)/x;
2015 : }
2016 : }
2017 :
2018 :
2019 0 : IPosition cursorShape(4, inx, 1, 1, 1);
2020 0 : IPosition axisPath(4, 0, 1, 2, 3);
2021 0 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
2022 0 : LatticeIterator<Complex> lix(*lattice, lsx);
2023 0 : for(lix.reset();!lix.atEnd();lix++) {
2024 0 : Int pol=lix.position()(2);
2025 0 : Int chan=lix.position()(3);
2026 0 : if(weights(pol, chan)!=0.0) {
2027 0 : Int iy=lix.position()(1);
2028 : //gridder->correctX1D(correction,iy);
2029 0 : for (Int ix=0;ix<nx;ix++) {
2030 0 : correction(ix)=1.0/(sincConvX(ix)*sincConvY(iy));
2031 : }
2032 0 : lix.rwVectorCursor()*=correction;
2033 0 : if(normalize) {
2034 0 : Complex rnorm(1.0/weights(pol,chan));
2035 0 : lix.rwCursor()*=rnorm;
2036 : }
2037 : }
2038 : else {
2039 0 : lix.woCursor()=0.0;
2040 : }
2041 : }
2042 0 : }
2043 :
2044 :
2045 :
2046 : //if(!isTiled)
2047 : {
2048 0 : LatticeLocker lock1 (*(image), FileLocker::Write);
2049 : // Check the section from the image BEFORE converting to a lattice
2050 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
2051 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
2052 0 : IPosition stride(4, 1);
2053 0 : IPosition trc(blc+image->shape()-stride);
2054 :
2055 : // Do the copy
2056 0 : IPosition start(4, 0);
2057 0 : image->put(griddedData(blc, trc));
2058 0 : }
2059 : }
2060 0 : if(!arrayLattice.null()) arrayLattice=0;
2061 0 : if(!lattice.null()) lattice=0;
2062 0 : griddedData.resize();
2063 0 : image->clearCache();
2064 0 : return *image;
2065 : }
2066 :
2067 : // Get weight image
2068 0 : void MosaicFT::getWeightImage(ImageInterface<Float>& weightImage,
2069 : Matrix<Float>& weights)
2070 : {
2071 :
2072 0 : logIO() << LogOrigin("MosaicFT", "getWeightImage") << LogIO::NORMAL;
2073 :
2074 0 : weights.resize(sumWeight.shape());
2075 0 : convertArray(weights,sumWeight);
2076 : /*
2077 : weightImage.copyData((LatticeExpr<Float>)
2078 : (iif((pbConvFunc_p->getFluxScaleImage()) > (0.0),
2079 : (*skyCoverage_p),0.0)));
2080 : */
2081 0 : weightImage.copyData(*skyCoverage_p);
2082 : //cerr << "getWeightIm " << max(sumWeight) << " " << max(skyCoverage_p->get()) << endl;
2083 :
2084 0 : skyCoverage_p->tempClose();
2085 :
2086 0 : }
2087 :
2088 0 : void MosaicFT::getFluxImage(ImageInterface<Float>& fluxImage) {
2089 :
2090 0 : if (stokes_p=="QU" || stokes_p=="IV"){
2091 0 : pbConvFunc_p->sliceFluxScale(2);
2092 : }
2093 0 : else if(stokes_p=="Q" || stokes_p=="U" || stokes_p=="V" || stokes_p=="I" ){
2094 0 : pbConvFunc_p->sliceFluxScale(1);
2095 : }
2096 0 : else if(stokes_p=="IQU"){
2097 0 : pbConvFunc_p->sliceFluxScale(3);
2098 : }
2099 0 : IPosition inShape=(pbConvFunc_p->getFluxScaleImage()).shape();
2100 0 : IPosition outShape=fluxImage.shape();
2101 0 : if(outShape==inShape){
2102 0 : fluxImage.copyData(pbConvFunc_p->getFluxScaleImage());
2103 : }
2104 0 : else if((outShape(0)==inShape(0)) && (outShape(1)==inShape(1))
2105 0 : && (outShape(2)==inShape(2))){
2106 : //case where CubeSkyEquation is chunking...copy the first pol-cube
2107 0 : IPosition cursorShape(4, inShape(0), inShape(1), inShape(2), 1);
2108 0 : IPosition axisPath(4, 0, 1, 2, 3);
2109 0 : LatticeStepper lsout(outShape, cursorShape, axisPath);
2110 0 : LatticeStepper lsin(inShape, cursorShape, axisPath);
2111 0 : LatticeIterator<Float> liout(fluxImage, lsout);
2112 0 : RO_LatticeIterator<Float> liin(pbConvFunc_p->getFluxScaleImage(), lsin);
2113 0 : liin.reset();
2114 0 : for(liout.reset();!liout.atEnd();liout++) {
2115 0 : if(inShape(2)==1)
2116 0 : liout.woMatrixCursor()=liin.matrixCursor();
2117 : else
2118 0 : liout.woCubeCursor()=liin.cubeCursor();
2119 : }
2120 :
2121 :
2122 0 : }
2123 : else{
2124 : //Should not reach here but the we're getting old
2125 0 : cout << "Bad case of shape mismatch in flux image shape" << endl;
2126 : }
2127 0 : }
2128 :
2129 0 : CountedPtr<TempImage<Float> >& MosaicFT::getConvWeightImage(){
2130 0 : if(!doneWeightImage_p)
2131 0 : finalizeToSky();
2132 0 : return skyCoverage_p;
2133 : }
2134 :
2135 6 : Bool MosaicFT::toRecord(String& error,
2136 : RecordInterface& outRec, Bool withImage, const String diskimage)
2137 : {
2138 : // Save the current MosaicFT object to an output state record
2139 6 : Bool retval = true;
2140 6 : if(!FTMachine::toRecord(error, outRec, withImage, diskimage))
2141 0 : return false;
2142 :
2143 6 : if(sj_p){
2144 6 : outRec.define("telescope", sj_p->telescope());
2145 : //cerr <<" Telescope " << sj_p->telescope() << endl;
2146 : }
2147 6 : outRec.define("uvscale", uvScale);
2148 6 : outRec.define("uvoffset", uvOffset);
2149 6 : outRec.define("cachesize", Int64(cachesize));
2150 6 : outRec.define("tilesize", tilesize);
2151 6 : outRec.define("maxabsdata", maxAbsData);
2152 6 : Vector<Int> center_loc(4), offset_loc(4);
2153 30 : for (Int k=0; k<4 ; k++){
2154 24 : center_loc(k)=centerLoc(k);
2155 24 : offset_loc(k)=offsetLoc(k);
2156 : }
2157 6 : outRec.define("centerloc", center_loc);
2158 6 : outRec.define("offsetloc", offset_loc);
2159 6 : outRec.define("usezero", usezero_p);
2160 6 : outRec.define("convfunc", convFunc);
2161 6 : outRec.define("weightconvfunc", weightConvFunc_p);
2162 6 : outRec.define("convsampling", convSampling);
2163 6 : outRec.define("convsize", convSize);
2164 6 : outRec.define("convsupport", convSupport);
2165 6 : outRec.define("convsupportplanes", convSupportPlanes_p);
2166 6 : outRec.define("convsizeplanes", convSizePlanes_p);
2167 6 : outRec.define("convRowMap", convRowMap_p);
2168 6 : outRec.define("stokes", stokes_p);
2169 6 : outRec.define("useconjconvfunc", useConjConvFunc_p);
2170 6 : outRec.define("usepointingtable", usePointingTable_p);
2171 6 : if(!pbConvFunc_p.null()){
2172 6 : Record subRec;
2173 : //cerr << "Doing pbconvrec " << endl;
2174 6 : pbConvFunc_p->toRecord(subRec);
2175 6 : outRec.defineRecord("pbconvfunc", subRec);
2176 6 : }
2177 :
2178 :
2179 6 : return retval;
2180 6 : }
2181 :
2182 0 : Bool MosaicFT::fromRecord(String& error,
2183 : const RecordInterface& inRec)
2184 : {
2185 0 : Bool retval = true;
2186 0 : pointingToImage=0;
2187 0 : doneWeightImage_p=false;
2188 0 : convWeightImage_p=nullptr;
2189 :
2190 0 : if(!FTMachine::fromRecord(error, inRec))
2191 0 : return false;
2192 0 : sj_p=0;
2193 0 : if(inRec.isDefined("telescope")){
2194 0 : String tel=inRec.asString("telescope");
2195 : PBMath::CommonPB pbtype;
2196 0 : Quantity freq(1e12, "Hz");// no useful band...just get default beam
2197 0 : String band="";
2198 0 : String pbname;
2199 0 : PBMath::whichCommonPBtoUse(tel, freq, band, pbtype, pbname);
2200 0 : if(pbtype != PBMath::UNKNOWN)
2201 0 : sj_p.reset(new VPSkyJones(tel,pbtype));
2202 0 : }
2203 :
2204 0 : inRec.get("name", machineName_p);
2205 0 : inRec.get("uvscale", uvScale);
2206 0 : inRec.get("uvoffset", uvOffset);
2207 0 : cachesize=inRec.asInt64("cachesize");
2208 0 : inRec.get("tilesize", tilesize);
2209 0 : inRec.get("maxabsdata", maxAbsData);
2210 0 : Vector<Int> center_loc(4), offset_loc(4);
2211 0 : inRec.get("centerloc", center_loc);
2212 0 : inRec.get("offsetloc", offset_loc);
2213 0 : uInt ndim4 = 4;
2214 0 : centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2),
2215 0 : center_loc(3));
2216 0 : offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2),
2217 0 : offset_loc(3));
2218 0 : imageCache=0; lattice=0; arrayLattice=0;
2219 0 : inRec.get("usezero", usezero_p);
2220 0 : inRec.get("convfunc", convFunc);
2221 0 : inRec.get("weightconvfunc", weightConvFunc_p);
2222 0 : inRec.get("convsampling", convSampling);
2223 0 : inRec.get("convsize", convSize);
2224 0 : inRec.get("convsupport", convSupport);
2225 0 : inRec.get("convsupportplanes", convSupportPlanes_p);
2226 0 : inRec.get("convsizeplanes", convSizePlanes_p);
2227 0 : inRec.get("convRowMap", convRowMap_p);
2228 0 : inRec.get("stokes", stokes_p);
2229 0 : inRec.get("useconjconvfunc", useConjConvFunc_p);
2230 0 : inRec.get("usepointingtable", usePointingTable_p);
2231 0 : if(inRec.isDefined("pbconvfunc")){
2232 0 : Record subRec=inRec.asRecord("pbconvfunc");
2233 0 : String elname=subRec.asString("name");
2234 : // if we are predicting only ...no need to estimate fluxscale
2235 0 : if(elname=="HetArrayConvFunc"){
2236 :
2237 0 : pbConvFunc_p=new HetArrayConvFunc(subRec, !toVis_p);
2238 : }
2239 : else{
2240 0 : pbConvFunc_p=new SimplePBConvFunc(subRec, !toVis_p);
2241 0 : if(!sj_p)
2242 0 : throw(AipsError("Failed to recovermosaic FTmachine;\n If you are seeing this message when try to get model vis \n then either try to reset the model or use scratch column for now"));
2243 : }
2244 0 : }
2245 : else{
2246 0 : pbConvFunc_p=0;
2247 : }
2248 0 : gridder.reset(nullptr);
2249 0 : return retval;
2250 0 : }
2251 :
2252 99486 : void MosaicFT::ok() {
2253 99486 : AlwaysAssert(image, AipsError);
2254 99486 : }
2255 :
2256 : // Make a plain straightforward honest-to-God image. This returns
2257 : // a complex image, without conversion to Stokes. The representation
2258 : // is that required for the visibilities.
2259 : //----------------------------------------------------------------------
2260 2 : void MosaicFT::makeImage(FTMachine::Type type,
2261 : vi::VisibilityIterator2& vi,
2262 : ImageInterface<Complex>& theImage,
2263 : Matrix<Float>& weight) {
2264 :
2265 :
2266 2 : logIO() << LogOrigin("MosaicFT", "makeImage") << LogIO::NORMAL;
2267 :
2268 2 : if(type==FTMachine::COVERAGE) {
2269 : logIO() << "Type COVERAGE not defined for Fourier transforms"
2270 0 : << LogIO::EXCEPTION;
2271 : }
2272 :
2273 :
2274 :
2275 : // Loop over all visibilities and pixels
2276 2 : vi::VisBuffer2* vb=vi.getVisBuffer();
2277 :
2278 : // Initialize put (i.e. transform to Sky) for this model
2279 2 : vi.origin();
2280 :
2281 2 : if(vb->polarizationFrame()==MSIter::Linear) {
2282 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
2283 : }
2284 : else {
2285 2 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
2286 : }
2287 :
2288 2 : initializeToSky(theImage,weight,*vb);
2289 : //This call is a NOP for all weighting schemes except for cube-briggs-perchanweightdensity
2290 2 : initBriggsWeightor(vi);
2291 : // Loop over the visibilities, putting VisBuffers
2292 10 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
2293 872 : for (vi.origin(); vi.more(); vi.next()) {
2294 :
2295 864 : switch(type) {
2296 0 : case FTMachine::RESIDUAL:
2297 0 : vb->setVisCube(vb->visCubeCorrected()-vb->visCubeModel());
2298 0 : put(*vb, -1, false);
2299 0 : break;
2300 0 : case FTMachine::MODEL:
2301 0 : vb->setVisCube(vb->visCubeModel());
2302 0 : put(*vb, -1, false);
2303 0 : break;
2304 0 : case FTMachine::CORRECTED:
2305 0 : vb->setVisCube(vb->visCubeCorrected());
2306 0 : put(*vb, -1, false);
2307 0 : break;
2308 864 : case FTMachine::PSF:
2309 864 : vb->setVisCube(Complex(1.0,0.0));
2310 864 : put(*vb, -1, true);
2311 864 : break;
2312 0 : case FTMachine::OBSERVED:
2313 : default:
2314 0 : put(*vb, -1, false);
2315 0 : break;
2316 : }
2317 : }
2318 : }
2319 2 : finalizeToSky();
2320 : // Normalize by dividing out weights, etc.
2321 2 : getImage(weight, true);
2322 2 : }
2323 :
2324 0 : Bool MosaicFT::getXYPos(const vi::VisBuffer2& vb, Int row) {
2325 :
2326 0 : MSColumns mscol(vb.ms());
2327 0 : const MSPointingColumns& act_mspc=mscol.pointing();
2328 0 : Int pointIndex=getIndex(act_mspc, vb.time()(row), vb.timeInterval()(row));
2329 0 : if((pointIndex<0)||pointIndex>=Int(act_mspc.time().nrow())) {
2330 : // ostringstream o;
2331 : // o << "Failed to find pointing information for time " <<
2332 : // MVTime(vb.time()(row)/86400.0);
2333 : // logIO_p << LogIO::DEBUGGING << String(o) << LogIO::POST;
2334 : // logIO_p << String(o) << LogIO::POST;
2335 :
2336 0 : worldPosMeas = mscol.field().phaseDirMeas(vb.fieldId()(0));
2337 : }
2338 : else {
2339 :
2340 0 : worldPosMeas=act_mspc.directionMeas(pointIndex);
2341 : // Make a machine to convert from the worldPosMeas to the output
2342 : // Direction Measure type for the relevant frame
2343 :
2344 :
2345 :
2346 : }
2347 :
2348 0 : if(!pointingToImage) {
2349 : // Set the frame - choose the first antenna. For e.g. VLBI, we
2350 : // will need to reset the frame per antenna
2351 0 : FTMachine::mLocation_p=mscol.antenna().positionMeas()(0);
2352 0 : mFrame_p=MeasFrame(MEpoch(Quantity(vb.time()(row), "s")),
2353 0 : FTMachine::mLocation_p);
2354 0 : MDirection::Ref outRef(directionCoord.directionType(), mFrame_p);
2355 0 : pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
2356 :
2357 0 : if(!pointingToImage) {
2358 0 : logIO_p << "Cannot make direction conversion machine" << LogIO::EXCEPTION;
2359 : }
2360 0 : }
2361 : else {
2362 0 : mFrame_p.resetEpoch(MEpoch(Quantity(vb.time()(row), "s")));
2363 : }
2364 :
2365 0 : worldPosMeas=(*pointingToImage)(worldPosMeas);
2366 :
2367 0 : Bool result=directionCoord.toPixel(xyPos, worldPosMeas);
2368 0 : if(!result) {
2369 0 : logIO_p << "Failed to find pixel location for "
2370 0 : << worldPosMeas.getAngle().getValue() << LogIO::EXCEPTION;
2371 0 : return false;
2372 : }
2373 0 : return result;
2374 :
2375 0 : }
2376 : // Get the index into the pointing table for this time. Note that the
2377 : // in the pointing table, TIME specifies the beginning of the spanned
2378 : // time range, whereas for the main table, TIME is the centroid.
2379 : // Note that the behavior for multiple matches is not specified! i.e.
2380 : // if there are multiple matches, the index returned depends on the
2381 : // history of previous matches. It is deterministic but not obvious.
2382 : // One could cure this by searching but it would be considerably
2383 : // costlier.
2384 0 : Int MosaicFT::getIndex(const MSPointingColumns& mspc, const Double& time,
2385 : const Double& /*interval*/) {
2386 0 : Int start=lastIndex_p;
2387 : // Search forwards
2388 0 : Int nrows=mspc.time().nrow();
2389 0 : if(nrows<1) {
2390 : // logIO_p << "No rows in POINTING table - cannot proceed" << LogIO::EXCEPTION;
2391 0 : return -1;
2392 : }
2393 0 : for (Int i=start;i<nrows;i++) {
2394 0 : Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
2395 : // If the interval in the pointing table is negative, use the last
2396 : // entry. Note that this may be invalid (-1) but in that case
2397 : // the calling routine will generate an error
2398 0 : if(mspc.interval()(i)<0.0) {
2399 0 : return lastIndex_p;
2400 : }
2401 : // Pointing table interval is specified so we have to do a match
2402 : else {
2403 : // Is the midpoint of this pointing table entry within the specified
2404 : // tolerance of the main table entry?
2405 0 : if(abs(midpoint-time) < (mspc.interval()(i)/2.0)) {
2406 0 : lastIndex_p=i;
2407 0 : return i;
2408 : }
2409 : }
2410 : }
2411 : // Search backwards
2412 0 : for (Int i=start;i>=0;i--) {
2413 0 : Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
2414 0 : if(mspc.interval()(i)<0.0) {
2415 0 : return lastIndex_p;
2416 : }
2417 : // Pointing table interval is specified so we have to do a match
2418 : else {
2419 : // Is the midpoint of this pointing table entry within the specified
2420 : // tolerance of the main table entry?
2421 0 : if(abs(midpoint-time) < (mspc.interval()(i)/2.0)) {
2422 0 : lastIndex_p=i;
2423 0 : return i;
2424 : }
2425 : }
2426 : }
2427 : // No match!
2428 0 : return -1;
2429 : }
2430 :
2431 :
2432 :
2433 :
2434 0 : void MosaicFT::addBeamCoverage(ImageInterface<Complex>& pbImage){
2435 :
2436 0 : CoordinateSystem cs(pbImage.coordinates());
2437 : // IPosition blc(4,0,0,0,0);
2438 : // IPosition trc(pbImage.shape());
2439 : // trc(0)=trc(0)-1;
2440 : // trc(1)=trc(1)-1;
2441 : // trc(2)=0;
2442 : // trc(3)=0;
2443 0 : WCBox *wbox= new WCBox(LCBox(pbImage.shape()), cs);
2444 0 : SubImage<Float> toAddTo(*skyCoverage_p, ImageRegion(wbox), true);
2445 0 : TempImage<Float> beamStokes(pbImage.shape(), cs);
2446 0 : StokesImageUtil::To(beamStokes, pbImage);
2447 : // toAddTo.copyData((LatticeExpr<Float>)(toAddTo + beamStokes ));
2448 0 : skyCoverage_p->copyData((LatticeExpr<Float>)(*skyCoverage_p + beamStokes ));
2449 :
2450 :
2451 0 : }
2452 :
2453 :
2454 :
2455 :
2456 :
2457 0 : String MosaicFT::name() const {
2458 0 : return machineName_p;
2459 : }
2460 :
2461 : } // REFIM ends
2462 : } //# NAMESPACE CASA - END
2463 :
2464 :
|
