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 207 : MosaicFT::MosaicFT(SkyJones* sj, MPosition mloc, String stokes,
94 : Long icachesize, Int itilesize,
95 207 : Bool usezero, Bool useDoublePrec, Bool useConjConvFunc, Bool usePointing)
96 207 : : FTMachine(), sj_p(sj),
97 207 : imageCache(0), cachesize(icachesize), tilesize(itilesize), gridder(nullptr),
98 207 : isTiled(false),
99 207 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
100 828 : mspc(0), msac(0), pointingToImage(0), usezero_p(usezero), convSampling(1),
101 621 : 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 207 : convSize=0;
104 207 : lastIndex_p=0;
105 207 : doneWeightImage_p=false;
106 207 : convWeightImage_p=0;
107 207 : pbConvFunc_p=new SimplePBConvFunc();
108 :
109 207 : mLocation_p=mloc;
110 207 : useDoubleGrid_p=useDoublePrec;
111 : // We should get rid of the ms dependence in the constructor
112 : // not used
113 207 : }
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 276 : MosaicFT& MosaicFT::operator=(const MosaicFT& other)
127 : {
128 276 : if(this!=&other) {
129 :
130 : //Do the base parameters
131 276 : FTMachine::operator=(other);
132 :
133 276 : convSampling=other.convSampling;
134 276 : sj_p=other.sj_p;
135 276 : imageCache=other.imageCache;
136 276 : cachesize=other.cachesize;
137 276 : tilesize=other.tilesize;
138 276 : isTiled=other.isTiled;
139 : //lattice=other.lattice;
140 276 : lattice=0;
141 : // arrayLattice=other.arrayLattice;
142 : // weightLattice=other.weightLattice;
143 : //if(arrayLattice) delete arrayLattice;
144 276 : arrayLattice=0;
145 : //if(weightLattice) delete weightLattice;
146 276 : weightLattice=0;
147 276 : maxAbsData=other.maxAbsData;
148 276 : centerLoc=other.centerLoc;
149 276 : offsetLoc=other.offsetLoc;
150 276 : pointingToImage=other.pointingToImage;
151 276 : usezero_p=other.usezero_p;
152 276 : doneWeightImage_p=other.doneWeightImage_p;
153 276 : pbConvFunc_p=other.pbConvFunc_p;
154 276 : stokes_p=other.stokes_p;
155 276 : if(!other.skyCoverage_p.null())
156 0 : skyCoverage_p=other.skyCoverage_p;
157 : else
158 276 : skyCoverage_p=0;
159 276 : if(other.convWeightImage_p !=0)
160 0 : convWeightImage_p=(TempImage<Complex> *)other.convWeightImage_p->cloneII();
161 : else
162 276 : convWeightImage_p=0;
163 276 : if(other.gridder==0)
164 276 : 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 276 : useConjConvFunc_p=other.useConjConvFunc_p;
174 276 : usePointingTable_p=other.usePointingTable_p;
175 276 : timemass_p=other.timemass_p;
176 276 : timegrid_p=other.timegrid_p;
177 276 : timedegrid_p=other.timedegrid_p;
178 : };
179 276 : return *this;
180 : };
181 :
182 : //----------------------------------------------------------------------
183 276 : MosaicFT::MosaicFT(const MosaicFT& other): FTMachine(),machineName_p("MosaicFT")
184 : {
185 276 : operator=(other);
186 276 : }
187 :
188 : //----------------------------------------------------------------------
189 : //void MosaicFT::setSharingFT(MosaicFT& otherFT){
190 : // otherFT_p=&otherFT;
191 : //}
192 370 : 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 370 : isTiled=false;
203 370 : nx = image->shape()(0);
204 370 : ny = image->shape()(1);
205 370 : npol = image->shape()(2);
206 370 : 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 370 : sumWeight.resize(npol, nchan);
217 :
218 370 : convSupport=0;
219 :
220 370 : uvScale.resize(2);
221 370 : uvScale=0.0;
222 370 : uvScale(0)=Float(nx)*image->coordinates().increment()(0);
223 370 : uvScale(1)=Float(ny)*image->coordinates().increment()(1);
224 :
225 370 : uvOffset.resize(2);
226 370 : uvOffset(0)=nx/2;
227 370 : uvOffset(1)=ny/2;
228 :
229 740 : gridder.reset(new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
230 370 : uvScale, uvOffset,
231 370 : "SF"));
232 :
233 : // Set up image cache needed for gridding.
234 370 : 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 370 : }
252 :
253 : // This is nasty, we should use CountedPointers here.
254 483 : MosaicFT::~MosaicFT() {
255 483 : if(imageCache) delete imageCache; imageCache=0;
256 : // if(arrayLattice) delete arrayLattice; arrayLattice=0;
257 483 : }
258 :
259 :
260 250 : void MosaicFT::setConvFunc(CountedPtr<SimplePBConvFunc>& pbconvFunc){
261 :
262 :
263 250 : pbConvFunc_p=pbconvFunc;
264 :
265 :
266 250 : }
267 :
268 69 : CountedPtr<SimplePBConvFunc>& MosaicFT::getConvFunc(){
269 69 : return pbConvFunc_p;
270 : }
271 :
272 85077 : void MosaicFT::findConvFunction(const ImageInterface<Complex>& iimage,
273 : const vi::VisBuffer2& vb, const Matrix<Double>& /*rotateduvw*/) {
274 :
275 :
276 : //oversample if image is small
277 : //But not more than 5000 pixels
278 85077 : convSampling=(max(nx, ny) < 50) ? 100: Int(ceil(5000.0/max(nx, ny)));
279 85077 : if(convSampling <1)
280 0 : convSampling=1;
281 85077 : if(pbConvFunc_p.null())
282 0 : pbConvFunc_p=new SimplePBConvFunc();
283 85077 : if(sj_p)
284 85077 : pbConvFunc_p->setSkyJones(sj_p.get());
285 : ////TEST for HetArray only for now
286 85077 : if(pbConvFunc_p->name()=="HetArrayConvFunc"){
287 71629 : if(convSampling <10)
288 30452 : convSampling=10;
289 71629 : AipsrcValue<Int>::find (convSampling, "mosaic.oversampling", 10);
290 : }
291 85077 : 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 85077 : pbConvFunc_p->findConvFunction(iimage, vb, convSampling, interpVisFreq_p, convFunc, weightConvFunc_p, convSizePlanes_p, convSupportPlanes_p,
299 85077 : 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 85077 : if(convSizePlanes_p.nelements() ==0)
304 0 : convSize=0;
305 : else
306 85077 : convSize=max(convSizePlanes_p);
307 85077 : if(convSupportPlanes_p.nelements() ==0)
308 0 : convSupport=0;
309 : else
310 85077 : convSupport=max(convSupportPlanes_p);
311 :
312 85077 : }
313 :
314 83 : void MosaicFT::initializeToVis(ImageInterface<Complex>& iimage,
315 : const vi::VisBuffer2& vb)
316 : {
317 83 : image=&iimage;
318 83 : toVis_p=true;
319 83 : ok();
320 :
321 : // if(convSize==0) {
322 83 : init(vb);
323 :
324 : // }
325 :
326 : // Initialize the maps for polarization and channel. These maps
327 : // translate visibility indices into image indices
328 83 : initMaps(vb);
329 83 : pbConvFunc_p->setVBUtil(vbutil_p);
330 83 : pbConvFunc_p->setUsePointing(usePointingTable_p);
331 : //make sure we rotate the first field too
332 83 : lastFieldId_p=-1;
333 83 : phaseShifter_p=new UVWMachine(*uvwMachine_p);
334 : //This is needed here as we need to know the grid correction before FFTing
335 83 : findConvFunction(*image, vb, vb.uvw());
336 :
337 83 : prepGridForDegrid();
338 :
339 83 : }
340 :
341 :
342 83 : void MosaicFT::prepGridForDegrid(){
343 :
344 : //For now isTiled=false
345 83 : isTiled=false;
346 83 : nx = image->shape()(0);
347 83 : ny = image->shape()(1);
348 83 : npol = image->shape()(2);
349 83 : nchan = image->shape()(3);
350 :
351 83 : IPosition gridShape(4, nx, ny, npol, nchan);
352 83 : griddedData.resize(gridShape);
353 83 : griddedData=Complex(0.0);
354 :
355 83 : IPosition stride(4, 1);
356 166 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
357 166 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
358 166 : IPosition trc(blc+image->shape()-stride);
359 :
360 83 : IPosition start(4, 0);
361 83 : griddedData(blc, trc) = image->getSlice(start, image->shape());
362 :
363 83 : image->clearCache();
364 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
365 83 : arrayLattice = new ArrayLattice<Complex>(griddedData);
366 83 : lattice=arrayLattice;
367 : {///UnDo the grid correction
368 83 : Int inx = lattice->shape()(0);
369 : //Int iny = lattice->shape()(1);
370 83 : Vector<Complex> correction(inx);
371 83 : correction=Complex(1.0, 0.0);
372 :
373 : // Int npixCorr= max(nx,ny);
374 83 : Vector<Float> sincConvX(nx);
375 67459 : for (Int ix=0;ix<nx;ix++) {
376 67376 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
377 67376 : if(ix==nx/2) {
378 83 : sincConvX(ix)=1.0;
379 : }
380 : else {
381 67293 : sincConvX(ix)=sin(x)/x;
382 : }
383 : }
384 83 : Vector<Float> sincConvY(ny);
385 67459 : for (Int ix=0;ix<ny;ix++) {
386 67376 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
387 67376 : if(ix==ny/2) {
388 83 : sincConvY(ix)=1.0;
389 : }
390 : else {
391 67293 : sincConvY(ix)=sin(x)/x;
392 : }
393 : }
394 :
395 83 : IPosition cursorShape(4, inx, 1, 1, 1);
396 83 : IPosition axisPath(4, 0, 1, 2, 3);
397 83 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
398 83 : LatticeIterator<Complex> lix(*lattice, lsx);
399 151915 : for(lix.reset();!lix.atEnd();lix++) {
400 : //Int pol=lix.position()(2);
401 : //Int chan=lix.position()(3);
402 :
403 151832 : Int iy=lix.position()(1);
404 : //gridder->correctX1D(correction,iy);
405 158599288 : for (Int ix=0;ix<nx;ix++) {
406 158447456 : correction(ix)=(sincConvX(ix)*sincConvY(iy));
407 : }
408 151832 : lix.rwVectorCursor()*=correction;
409 :
410 : }
411 83 : }
412 :
413 :
414 83 : logIO() << LogIO::DEBUGGING << "Starting FFT of image" << LogIO::POST;
415 : // Now do the FFT2D in place
416 83 : 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 83 : logIO() << LogIO::DEBUGGING << "Finished FFT" << LogIO::POST;
432 :
433 :
434 :
435 83 : }
436 :
437 :
438 83 : void MosaicFT::finalizeToVis()
439 : {
440 83 : logIO() << LogOrigin("MosaicFT", "finalizeToVis") << LogIO::NORMAL;
441 83 : logIO()<< LogIO::NORMAL2 << "Time degrid " << timedegrid_p << LogIO::POST;
442 83 : timedegrid_p=0.0;
443 :
444 83 : if(!arrayLattice.null()) arrayLattice=0;
445 83 : if(!lattice.null()) lattice=0;
446 83 : 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 83 : if(pointingToImage)
462 0 : delete pointingToImage;
463 83 : pointingToImage=0;
464 83 : }
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 287 : void MosaicFT::initializeToSky(ImageInterface<Complex>& iimage,
475 : Matrix<Float>& weight,
476 : const vi::VisBuffer2& vb)
477 : {
478 : // image always points to the image
479 287 : image=&iimage;
480 287 : toVis_p=False;
481 : // if(convSize==0) {
482 287 : init(vb);
483 :
484 : // }
485 :
486 : // Initialize the maps for polarization and channel. These maps
487 : // translate visibility indices into image indices
488 287 : initMaps(vb);
489 287 : pbConvFunc_p->setVBUtil(vbutil_p);
490 287 : pbConvFunc_p->setUsePointing(usePointingTable_p);
491 : //make sure we rotate the first field too
492 287 : lastFieldId_p=-1;
493 287 : 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 287 : isTiled=false;
504 287 : nx = image->shape()(0);
505 287 : ny = image->shape()(1);
506 287 : npol = image->shape()(2);
507 287 : nchan = image->shape()(3);
508 :
509 287 : sumWeight=0.0;
510 287 : weight.resize(sumWeight.shape());
511 287 : weight=0.0;
512 :
513 287 : 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 287 : IPosition gridShape(4, nx, ny, npol, nchan);
537 287 : if(!useDoubleGrid_p) {
538 0 : griddedData.resize(gridShape);
539 0 : griddedData=Complex(0.0);
540 : }
541 : else {
542 287 : griddedData2.resize(gridShape);
543 287 : griddedData2=DComplex(0.0);
544 : }
545 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
546 : //arrayLattice = new ArrayLattice<Complex>(griddedData);
547 : //lattice=arrayLattice;
548 :
549 287 : if( !doneWeightImage_p && (convWeightImage_p==0)){
550 :
551 :
552 :
553 408 : convWeightImage_p=new TempImage<Complex> (iimage.shape(),
554 204 : iimage.coordinates());
555 204 : 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 204 : if(useDoubleGrid_p){
564 204 : griddedWeight2.resize(gridShape);
565 204 : griddedWeight2=DComplex(0.0);
566 : }
567 : else{
568 0 : griddedWeight=Complex(0.0);
569 : }
570 : //if(weightLattice) delete weightLattice; weightLattice=0;
571 204 : weightLattice = new ArrayLattice<Complex>(griddedWeight);
572 :
573 : }
574 :
575 287 : }
576 :
577 : //cerr << "initializetosky lastfield " << lastFieldId_p << endl;
578 : // AlwaysAssert(lattice, AipsError);
579 :
580 287 : }
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 287 : void MosaicFT::finalizeToSky()
591 : {
592 287 : logIO() << LogOrigin("MosaicFT", "finalizeToSky") << LogIO::NORMAL;
593 287 : logIO() << LogIO::NORMAL2 << "time to massage data " << timemass_p << LogIO::POST;
594 287 : logIO() << LogIO::NORMAL2<< "time gridding " << timegrid_p << LogIO::POST;
595 287 : timemass_p=0.0;
596 287 : 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 287 : if(!doneWeightImage_p){
613 204 : if(useDoubleGrid_p){
614 204 : convertArray(griddedWeight, griddedWeight2);
615 : //Don't need the double-prec grid anymore...
616 204 : griddedWeight2.resize();
617 : }
618 204 : ft_p.c2cFFT(*weightLattice, false);
619 : //Get the stokes right
620 204 : CoordinateSystem coords=convWeightImage_p->coordinates();
621 204 : Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
622 204 : Int npol=1;
623 204 : Vector<Int> whichStokes(npol);
624 204 : if(stokes_p=="I" || stokes_p=="RR" || stokes_p=="LL" ||stokes_p=="XX"
625 204 : || stokes_p=="YY" || stokes_p=="Q" || stokes_p=="U" || stokes_p=="V"){
626 204 : npol=1;
627 204 : 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 0 : else if(stokes_p=="IV"){
633 0 : npol=2;
634 0 : whichStokes.resize(2);
635 0 : whichStokes(0)=Stokes::I;
636 0 : whichStokes(1)=Stokes::V;
637 : }
638 0 : else if(stokes_p=="QU"){
639 0 : npol=2;
640 0 : whichStokes.resize(2);
641 0 : whichStokes(0)=Stokes::Q;
642 0 : whichStokes(1)=Stokes::U;
643 : }
644 0 : 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 0 : 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 0 : 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 0 : else if(stokes_p=="IQUV"){
664 0 : npol=4;
665 0 : whichStokes.resize(4);
666 0 : whichStokes(0)=Stokes::I;
667 0 : whichStokes(1)=Stokes::Q;
668 0 : whichStokes(2)=Stokes::U;
669 0 : whichStokes(3)=Stokes::V;
670 : }
671 :
672 204 : StokesCoordinate newStokesCoord(whichStokes);
673 204 : coords.replaceCoordinate(newStokesCoord, stokesIndex);
674 204 : IPosition imshp=convWeightImage_p->shape();
675 204 : imshp(2)=npol;
676 :
677 :
678 204 : skyCoverage_p=new TempImage<Float> (imshp, coords,1.0);
679 408 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
680 408 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
681 204 : IPosition stride(4, 1);
682 408 : IPosition trc(blc+image->shape()-stride);
683 :
684 : // Do the copy
685 204 : IPosition start(4, 0);
686 204 : convWeightImage_p->put(griddedWeight(blc, trc));
687 204 : StokesImageUtil::ToStokesPSF(*skyCoverage_p, *convWeightImage_p);
688 204 : 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 0 : blc(0)=0; blc(1)=0; blc(3)=0;blc(2)=0;
692 0 : trc=skyCoverage_p->shape()-stride;
693 0 : trc(2)=0;
694 0 : SubImage<Float> isubim(*skyCoverage_p, Slicer(blc, trc, Slicer::endIsLast));
695 0 : for (Int k=1; k < npol; ++k){
696 0 : blc(2)=k; trc(2)=k;
697 0 : SubImage<Float> quvsubim(*skyCoverage_p, Slicer(blc, trc, Slicer::endIsLast), true);
698 0 : quvsubim.copyData(isubim);
699 0 : }
700 :
701 0 : }
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 204 : pbConvFunc_p->setWeightImage(skyCoverage_p);
706 204 : if(convWeightImage_p) delete convWeightImage_p;
707 204 : convWeightImage_p=0;
708 204 : 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 204 : }
719 :
720 287 : if(!weightLattice.null()) weightLattice=0;
721 287 : griddedWeight.resize();
722 287 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
723 287 : }
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 65161 : void MosaicFT::put(const vi::VisBuffer2& vb, Int row, Bool dopsf,
994 : FTMachine::Type type)
995 : {
996 :
997 :
998 :
999 :
1000 65161 : Timer tim;
1001 65161 : tim.mark();
1002 :
1003 65161 : matchChannel(vb);
1004 :
1005 :
1006 : //cerr << "CHANMAP " << chanMap << endl;
1007 : //No point in reading data if its not matching in frequency
1008 65161 : if(max(chanMap)==-1)
1009 0 : return;
1010 :
1011 : //const Matrix<Float> *imagingweight;
1012 : //imagingweight=&(vb.imagingWeight());
1013 65161 : Matrix<Float> imagingweight;
1014 65161 : getImagingWeight(imagingweight, vb);
1015 :
1016 65161 : if(dopsf) type=FTMachine::PSF;
1017 :
1018 65161 : Cube<Complex> data;
1019 : //Fortran gridder need the flag as ints
1020 65161 : Cube<Int> flags;
1021 65161 : Matrix<Float> elWeight;
1022 65161 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
1023 :
1024 :
1025 :
1026 : Bool iswgtCopy;
1027 : const Float *wgtStorage;
1028 65161 : wgtStorage=elWeight.getStorage(iswgtCopy);
1029 :
1030 :
1031 :
1032 :
1033 : Bool isCopy;
1034 65161 : const Complex *datStorage=0;
1035 :
1036 : // cerr << "dopsf " << dopsf << " isWeightCopy " << iswgtCopy << " " << wgtStorage<< endl;
1037 65161 : if(!dopsf)
1038 39896 : datStorage=data.getStorage(isCopy);
1039 :
1040 :
1041 : // If row is -1 then we pass through all rows
1042 : Int startRow, endRow, nRow;
1043 65161 : if (row==-1) {
1044 65161 : nRow=vb.nRows();
1045 65161 : startRow=0;
1046 65161 : 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 65161 : Matrix<Double> uvw(negateUV(vb));
1056 65161 : Vector<Double> dphase(vb.nRows());
1057 65161 : dphase=0.0;
1058 :
1059 65161 : doUVWRotation_p=true;
1060 65161 : girarUVW(uvw, dphase, vb);
1061 65161 : 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 65161 : findConvFunction(*image, vb, uvw);
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 65161 : 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 65161 : Int idopsf=0;
1086 65161 : if(dopsf) idopsf=1;
1087 :
1088 :
1089 65161 : Vector<Int> rowFlags(vb.nRows());
1090 65161 : rowFlags=0;
1091 65161 : rowFlags(vb.flagRow())=true;
1092 65161 : if(!usezero_p) {
1093 15088945 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1094 15023784 : 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 65161 : const IPosition& fs=flags.shape();
1126 : //cerr << "flags shape " << fs << endl;
1127 65161 : std::vector<Int>s(fs.begin(), fs.end());
1128 65161 : Int nvp=s[0];
1129 65161 : Int nvc=s[1];
1130 65161 : Int nvisrow=s[2];
1131 65161 : Int csamp=convSampling;
1132 : Bool uvwcopy;
1133 65161 : const Double *uvwstor=uvw.getStorage(uvwcopy);
1134 : Bool gridcopy;
1135 : Bool convcopy;
1136 : Bool wconvcopy;
1137 65161 : const Complex *convstor=convFunc.getStorage(convcopy);
1138 65161 : const Complex *wconvstor=weightConvFunc_p.getStorage(wconvcopy);
1139 65161 : Int nPolConv=convFunc.shape()[2];
1140 65161 : Int nChanConv=convFunc.shape()[3];
1141 65161 : Int nConvFunc=convFunc.shape()(4);
1142 : Bool weightcopy;
1143 : ////////**************************
1144 65161 : Cube<Int> loc(2, nvc, nRow);
1145 65161 : Cube<Int> off(2, nvc, nRow);
1146 65161 : Matrix<Complex> phasor(nvc, nRow);
1147 : Bool delphase;
1148 65161 : Complex * phasorstor=phasor.getStorage(delphase);
1149 65161 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1150 65161 : const Double * scalestor=uvScale.getStorage(del);
1151 65161 : const Double * offsetstor=uvOffset.getStorage(del);
1152 65161 : Int * locstor=loc.getStorage(del);
1153 65161 : Int * offstor=off.getStorage(del);
1154 65161 : const Double *dpstor=dphase.getStorage(del);
1155 : Int irow;
1156 65161 : Int nth=1;
1157 : #ifdef _OPENMP
1158 65161 : if(numthreads_p >0){
1159 0 : nth=min(numthreads_p, omp_get_max_threads());
1160 : }
1161 : else{
1162 65161 : nth= omp_get_max_threads();
1163 : }
1164 : //nth=min(4,nth);
1165 : #endif
1166 65161 : Double cinv=Double(1.0)/C::c;
1167 :
1168 65161 : Int dow=0;
1169 65161 : #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 65161 : timemass_p +=tim.real();
1184 : Int ixsub, iysub, icounter;
1185 65161 : ixsub=1;
1186 65161 : iysub=1;
1187 : //////***********************DEBUGGING
1188 : //nth=1;
1189 : ////////***************
1190 65161 : if (nth >3){
1191 65161 : ixsub=8;
1192 65161 : iysub=8;
1193 : }
1194 0 : else if(nth >1){
1195 0 : ixsub=2;
1196 0 : iysub=2;
1197 : }
1198 65161 : Int rbeg=startRow+1;
1199 65161 : Int rend=endRow+1;
1200 65161 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
1201 65161 : Vector<Double *> swgtptr(ixsub*iysub);
1202 65161 : Vector<Bool> swgtdel(ixsub*iysub);
1203 4235465 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1204 4170304 : sumwgt[icounter].resize(sumWeight.shape());
1205 4170304 : sumwgt[icounter].set(0.0);
1206 4170304 : swgtptr[icounter]=sumwgt[icounter].getStorage(swgtdel(icounter));
1207 : }
1208 : //cerr << "done thread " << doneThreadPartition_p << " " << ixsub*iysub << endl;
1209 65161 : if(doneThreadPartition_p < 0){
1210 204 : xsect_p.resize(ixsub*iysub);
1211 204 : ysect_p.resize(ixsub*iysub);
1212 204 : nxsect_p.resize(ixsub*iysub);
1213 204 : nysect_p.resize(ixsub*iysub);
1214 13260 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1215 13056 : findGridSector(nx, ny, ixsub, iysub, 0, 0, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
1216 : }
1217 : }
1218 65161 : Vector<Int> xsect, ysect, nxsect, nysect;
1219 65161 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
1220 : //cerr << xsect.shape() << " " << xsect << endl;
1221 65161 : const Int* pmapstor=polMap.getStorage(del);
1222 65161 : const Int* cmapstor=chanMap.getStorage(del);
1223 : // Dummy sumwt for gridweight part
1224 65161 : Matrix<Double> dumSumWeight(npol, nchan);
1225 65161 : dumSumWeight=sumWeight;
1226 : Bool isDSWC;
1227 65161 : Double *dsumwtstor=dumSumWeight.getStorage(isDSWC);
1228 65161 : Int nc=nchan;
1229 65161 : Int np=npol;
1230 65161 : Int nxp=nx;
1231 65161 : Int nyp=ny;
1232 65161 : Int csize=convSize;
1233 65161 : const Int * flagstor=flags.getStorage(del);
1234 65161 : const Int * rowflagstor=rowFlags.getStorage(del);
1235 65161 : Int csupp=convSupport;
1236 65161 : const Int *convrowmapstor=convRowMap_p.getStorage(del);
1237 65161 : const Int *convchanmapstor=convChanMap_p.getStorage(del);
1238 65161 : const Int *convpolmapstor=convPolMap_p.getStorage(del);
1239 : ///
1240 :
1241 :
1242 : ////////***************************
1243 65161 : tim.mark();
1244 :
1245 65161 : if(useDoubleGrid_p) {
1246 65161 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
1247 :
1248 65161 : #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 4235465 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1302 4170304 : sumwgt[icounter].putStorage(swgtptr[icounter],swgtdel[icounter]);
1303 4170304 : sumWeight=sumWeight+sumwgt[icounter];
1304 : }
1305 :
1306 : //cerr << "SUMWEIG " << sumWeight << endl;
1307 65161 : griddedData2.putStorage(gridstor, gridcopy);
1308 65161 : if(dopsf && (nth >4))
1309 10440 : tweakGridSector(nx, ny, ixsub, iysub);
1310 65161 : timegrid_p+=tim.real();
1311 65161 : tim.mark();
1312 65161 : 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 43659 : DComplex *gridwgtstor=griddedWeight2.getStorage(weightcopy);
1316 43659 : 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 43659 : griddedWeight2.putStorage(gridwgtstor, weightcopy);
1324 :
1325 : }
1326 65161 : 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 65161 : convFunc.freeStorage(convstor, convcopy);
1409 65161 : weightConvFunc_p.freeStorage(wconvstor, wconvcopy);
1410 65161 : dumSumWeight.putStorage(dsumwtstor, isDSWC);
1411 : //cerr << "dumSumwe " << dumSumWeight << endl;
1412 65161 : uvw.freeStorage(uvwstor, uvwcopy);
1413 65161 : if(!dopsf)
1414 39896 : data.freeStorage(datStorage, isCopy);
1415 :
1416 65161 : elWeight.freeStorage(wgtStorage,iswgtCopy);
1417 :
1418 :
1419 :
1420 :
1421 65161 : }
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 19833 : 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 19833 : if (row==-1) {
1618 19833 : nRow=vb.nRows();
1619 19833 : startRow=0;
1620 19833 : 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 19833 : matchChannel(vb);
1633 :
1634 : //No point in reading data if its not matching in frequency
1635 19833 : if(max(chanMap)==-1)
1636 0 : return;
1637 :
1638 : // Get the uvws in a form that Fortran can use
1639 19833 : Matrix<Double> uvw(negateUV(vb));
1640 19833 : Vector<Double> dphase(vb.nRows());
1641 19833 : dphase=0.0;
1642 :
1643 19833 : doUVWRotation_p=true;
1644 19833 : girarUVW(uvw, dphase, vb);
1645 19833 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1646 :
1647 :
1648 :
1649 :
1650 19833 : Cube<Complex> data;
1651 19833 : Cube<Int> flags;
1652 19833 : 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 19833 : findConvFunction(*image, vb, uvw);
1659 :
1660 : // no valid pointing in this buffer
1661 19833 : if(convSupport <= 0)
1662 0 : return;
1663 : Complex *datStorage;
1664 : Bool isCopy;
1665 19833 : datStorage=data.getStorage(isCopy);
1666 :
1667 :
1668 19833 : Vector<Int> rowFlags(vb.nRows());
1669 19833 : rowFlags=0;
1670 19833 : rowFlags(vb.flagRow())=true;
1671 19833 : if(!usezero_p) {
1672 2981715 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1673 2961882 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
1674 : }
1675 : }
1676 19833 : Int nvp=data.shape()[0];
1677 19833 : Int nvc=data.shape()[1];
1678 19833 : Int nvisrow=data.shape()[2];
1679 19833 : Int csamp=convSampling;
1680 19833 : Int csize=convSize;
1681 19833 : Int csupp=convSupport;
1682 19833 : Int nc=nchan;
1683 19833 : Int np=npol;
1684 19833 : Int nxp=nx;
1685 19833 : Int nyp=ny;
1686 : Bool uvwcopy;
1687 19833 : const Double *uvwstor=uvw.getStorage(uvwcopy);
1688 19833 : Int nPolConv=convFunc.shape()[2];
1689 19833 : Int nChanConv=convFunc.shape()[3];
1690 19833 : Int nConvFunc=convFunc.shape()(4);
1691 : ////////**************************
1692 19833 : Cube<Int> loc(2, nvc, nRow);
1693 19833 : Cube<Int> off(2, nvc, nRow);
1694 19833 : Matrix<Complex> phasor(nvc, nRow);
1695 : Bool delphase;
1696 : Bool del;
1697 19833 : const Int* pmapstor=polMap.getStorage(del);
1698 19833 : const Int* cmapstor=chanMap.getStorage(del);
1699 19833 : Complex * phasorstor=phasor.getStorage(delphase);
1700 19833 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1701 19833 : const Double * scalestor=uvScale.getStorage(del);
1702 19833 : const Double * offsetstor=uvOffset.getStorage(del);
1703 19833 : const Int * flagstor=flags.getStorage(del);
1704 19833 : const Int * rowflagstor=rowFlags.getStorage(del);
1705 19833 : Int * locstor=loc.getStorage(del);
1706 19833 : Int * offstor=off.getStorage(del);
1707 19833 : const Double *dpstor=dphase.getStorage(del);
1708 19833 : const Int *convrowmapstor=convRowMap_p.getStorage(del);
1709 19833 : const Int *convchanmapstor=convChanMap_p.getStorage(del);
1710 19833 : const Int *convpolmapstor=convPolMap_p.getStorage(del);
1711 : ////////***************************
1712 :
1713 : Int irow;
1714 19833 : Int nth=1;
1715 : #ifdef _OPENMP
1716 19833 : if(numthreads_p >0){
1717 0 : nth=min(numthreads_p, omp_get_max_threads());
1718 : }
1719 : else{
1720 19833 : nth= omp_get_max_threads();
1721 : }
1722 : //nth=min(4,nth);
1723 : #endif
1724 :
1725 19833 : Timer tim;
1726 19833 : tim.mark();
1727 :
1728 19833 : Int dow=0;
1729 19833 : Double cinv=Double(1.0)/C::c;
1730 19833 : #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 19833 : Int rbeg=startRow+1;
1744 19833 : Int rend=endRow+1;
1745 19833 : Int npart=nth;
1746 :
1747 : Bool gridcopy;
1748 19833 : const Complex *gridstor=griddedData.getStorage(gridcopy);
1749 : Bool convcopy;
1750 : ////Degridding needs the conjugate ...doing it here
1751 19833 : Array<Complex> conjConvFunc=conj(convFunc);
1752 19833 : const Complex *convstor=conjConvFunc.getStorage(convcopy);
1753 19833 : Int ix=0;
1754 19833 : #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 19833 : data.putStorage(datStorage, isCopy);
1792 19833 : griddedData.freeStorage(gridstor, gridcopy);
1793 19833 : convFunc.freeStorage(convstor, convcopy);
1794 :
1795 19833 : timedegrid_p+=tim.real();
1796 :
1797 19833 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1798 19833 : }
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 85077 : void MosaicFT::ok() {
2253 85077 : AlwaysAssert(image, AipsError);
2254 85077 : }
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 :
|
