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 208 : MosaicFT::MosaicFT(SkyJones* sj, MPosition mloc, String stokes,
94 : Long icachesize, Int itilesize,
95 208 : Bool usezero, Bool useDoublePrec, Bool useConjConvFunc, Bool usePointing)
96 208 : : FTMachine(), sj_p(sj),
97 208 : imageCache(0), cachesize(icachesize), tilesize(itilesize), gridder(nullptr),
98 208 : isTiled(false),
99 208 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
100 832 : mspc(0), msac(0), pointingToImage(0), usezero_p(usezero), convSampling(1),
101 624 : 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 208 : convSize=0;
104 208 : lastIndex_p=0;
105 208 : doneWeightImage_p=false;
106 208 : convWeightImage_p=0;
107 208 : pbConvFunc_p=new SimplePBConvFunc();
108 :
109 208 : mLocation_p=mloc;
110 208 : useDoubleGrid_p=useDoublePrec;
111 : // We should get rid of the ms dependence in the constructor
112 : // not used
113 208 : }
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 277 : MosaicFT& MosaicFT::operator=(const MosaicFT& other)
127 : {
128 277 : if(this!=&other) {
129 :
130 : //Do the base parameters
131 277 : FTMachine::operator=(other);
132 :
133 277 : convSampling=other.convSampling;
134 277 : sj_p=other.sj_p;
135 277 : imageCache=other.imageCache;
136 277 : cachesize=other.cachesize;
137 277 : tilesize=other.tilesize;
138 277 : isTiled=other.isTiled;
139 : //lattice=other.lattice;
140 277 : lattice=0;
141 : // arrayLattice=other.arrayLattice;
142 : // weightLattice=other.weightLattice;
143 : //if(arrayLattice) delete arrayLattice;
144 277 : arrayLattice=0;
145 : //if(weightLattice) delete weightLattice;
146 277 : weightLattice=0;
147 277 : maxAbsData=other.maxAbsData;
148 277 : centerLoc=other.centerLoc;
149 277 : offsetLoc=other.offsetLoc;
150 277 : pointingToImage=other.pointingToImage;
151 277 : usezero_p=other.usezero_p;
152 277 : doneWeightImage_p=other.doneWeightImage_p;
153 277 : pbConvFunc_p=other.pbConvFunc_p;
154 277 : stokes_p=other.stokes_p;
155 277 : if(!other.skyCoverage_p.null())
156 0 : skyCoverage_p=other.skyCoverage_p;
157 : else
158 277 : skyCoverage_p=0;
159 277 : if(other.convWeightImage_p !=0)
160 0 : convWeightImage_p=(TempImage<Complex> *)other.convWeightImage_p->cloneII();
161 : else
162 277 : convWeightImage_p=0;
163 277 : if(other.gridder==0)
164 277 : 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 277 : useConjConvFunc_p=other.useConjConvFunc_p;
174 277 : usePointingTable_p=other.usePointingTable_p;
175 277 : timemass_p=other.timemass_p;
176 277 : timegrid_p=other.timegrid_p;
177 277 : timedegrid_p=other.timedegrid_p;
178 : };
179 277 : return *this;
180 : };
181 :
182 : //----------------------------------------------------------------------
183 277 : MosaicFT::MosaicFT(const MosaicFT& other): FTMachine(),machineName_p("MosaicFT")
184 : {
185 277 : operator=(other);
186 277 : }
187 :
188 : //----------------------------------------------------------------------
189 : //void MosaicFT::setSharingFT(MosaicFT& otherFT){
190 : // otherFT_p=&otherFT;
191 : //}
192 372 : 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 372 : isTiled=false;
203 372 : nx = image->shape()(0);
204 372 : ny = image->shape()(1);
205 372 : npol = image->shape()(2);
206 372 : 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 372 : sumWeight.resize(npol, nchan);
217 :
218 372 : convSupport=0;
219 :
220 372 : uvScale.resize(2);
221 372 : uvScale=0.0;
222 372 : uvScale(0)=Float(nx)*image->coordinates().increment()(0);
223 372 : uvScale(1)=Float(ny)*image->coordinates().increment()(1);
224 :
225 372 : uvOffset.resize(2);
226 372 : uvOffset(0)=nx/2;
227 372 : uvOffset(1)=ny/2;
228 :
229 744 : gridder.reset(new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
230 372 : uvScale, uvOffset,
231 372 : "SF"));
232 :
233 : // Set up image cache needed for gridding.
234 372 : 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 372 : }
252 :
253 : // This is nasty, we should use CountedPointers here.
254 485 : MosaicFT::~MosaicFT() {
255 485 : if(imageCache) delete imageCache; imageCache=0;
256 : // if(arrayLattice) delete arrayLattice; arrayLattice=0;
257 485 : }
258 :
259 :
260 251 : void MosaicFT::setConvFunc(CountedPtr<SimplePBConvFunc>& pbconvFunc){
261 :
262 :
263 251 : pbConvFunc_p=pbconvFunc;
264 :
265 :
266 251 : }
267 :
268 69 : CountedPtr<SimplePBConvFunc>& MosaicFT::getConvFunc(){
269 69 : return pbConvFunc_p;
270 : }
271 :
272 86517 : 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 86517 : convSampling=(max(nx, ny) < 50) ? 100: Int(ceil(5000.0/max(nx, ny)));
279 86517 : if(convSampling <1)
280 0 : convSampling=1;
281 86517 : if(pbConvFunc_p.null())
282 0 : pbConvFunc_p=new SimplePBConvFunc();
283 86517 : if(sj_p)
284 86517 : pbConvFunc_p->setSkyJones(sj_p.get());
285 : ////TEST for HetArray only for now
286 86517 : if(pbConvFunc_p->name()=="HetArrayConvFunc"){
287 73069 : if(convSampling <10)
288 30452 : convSampling=10;
289 73069 : AipsrcValue<Int>::find (convSampling, "mosaic.oversampling", 10);
290 : }
291 86517 : 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 86517 : pbConvFunc_p->findConvFunction(iimage, vb, convSampling, interpVisFreq_p, convFunc, weightConvFunc_p, convSizePlanes_p, convSupportPlanes_p,
299 86517 : 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 86517 : if(convSizePlanes_p.nelements() ==0)
304 0 : convSize=0;
305 : else
306 86517 : convSize=max(convSizePlanes_p);
307 86517 : if(convSupportPlanes_p.nelements() ==0)
308 0 : convSupport=0;
309 : else
310 86517 : convSupport=max(convSupportPlanes_p);
311 :
312 86517 : }
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 289 : void MosaicFT::initializeToSky(ImageInterface<Complex>& iimage,
475 : Matrix<Float>& weight,
476 : const vi::VisBuffer2& vb)
477 : {
478 : // image always points to the image
479 289 : image=&iimage;
480 289 : toVis_p=False;
481 : // if(convSize==0) {
482 289 : init(vb);
483 :
484 : // }
485 :
486 : // Initialize the maps for polarization and channel. These maps
487 : // translate visibility indices into image indices
488 289 : initMaps(vb);
489 289 : pbConvFunc_p->setVBUtil(vbutil_p);
490 289 : pbConvFunc_p->setUsePointing(usePointingTable_p);
491 : //make sure we rotate the first field too
492 289 : lastFieldId_p=-1;
493 289 : 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 289 : isTiled=false;
504 289 : nx = image->shape()(0);
505 289 : ny = image->shape()(1);
506 289 : npol = image->shape()(2);
507 289 : nchan = image->shape()(3);
508 :
509 289 : sumWeight=0.0;
510 289 : weight.resize(sumWeight.shape());
511 289 : weight=0.0;
512 :
513 289 : 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 289 : IPosition gridShape(4, nx, ny, npol, nchan);
537 289 : if(!useDoubleGrid_p) {
538 0 : griddedData.resize(gridShape);
539 0 : griddedData=Complex(0.0);
540 : }
541 : else {
542 289 : griddedData2.resize(gridShape);
543 289 : griddedData2=DComplex(0.0);
544 : }
545 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
546 : //arrayLattice = new ArrayLattice<Complex>(griddedData);
547 : //lattice=arrayLattice;
548 :
549 289 : if( !doneWeightImage_p && (convWeightImage_p==0)){
550 :
551 :
552 :
553 410 : convWeightImage_p=new TempImage<Complex> (iimage.shape(),
554 205 : iimage.coordinates());
555 205 : 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 205 : if(useDoubleGrid_p){
564 205 : griddedWeight2.resize(gridShape);
565 205 : griddedWeight2=DComplex(0.0);
566 : }
567 : else{
568 0 : griddedWeight=Complex(0.0);
569 : }
570 : //if(weightLattice) delete weightLattice; weightLattice=0;
571 205 : weightLattice = new ArrayLattice<Complex>(griddedWeight);
572 :
573 : }
574 :
575 289 : }
576 :
577 : //cerr << "initializetosky lastfield " << lastFieldId_p << endl;
578 : // AlwaysAssert(lattice, AipsError);
579 :
580 289 : }
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 289 : void MosaicFT::finalizeToSky()
591 : {
592 289 : logIO() << LogOrigin("MosaicFT", "finalizeToSky") << LogIO::NORMAL;
593 289 : logIO() << LogIO::NORMAL2 << "time to massage data " << timemass_p << LogIO::POST;
594 289 : logIO() << LogIO::NORMAL2<< "time gridding " << timegrid_p << LogIO::POST;
595 289 : timemass_p=0.0;
596 289 : 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 289 : if(!doneWeightImage_p){
613 205 : if(useDoubleGrid_p){
614 205 : convertArray(griddedWeight, griddedWeight2);
615 : //Don't need the double-prec grid anymore...
616 205 : griddedWeight2.resize();
617 : }
618 205 : ft_p.c2cFFT(*weightLattice, false);
619 : //Get the stokes right
620 205 : CoordinateSystem coords=convWeightImage_p->coordinates();
621 205 : Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
622 205 : Int npol=1;
623 205 : Vector<Int> whichStokes(npol);
624 205 : if(stokes_p=="I" || stokes_p=="RR" || stokes_p=="LL" ||stokes_p=="XX"
625 205 : || stokes_p=="YY" || stokes_p=="Q" || stokes_p=="U" || stokes_p=="V"){
626 205 : npol=1;
627 205 : 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 205 : StokesCoordinate newStokesCoord(whichStokes);
673 205 : coords.replaceCoordinate(newStokesCoord, stokesIndex);
674 205 : IPosition imshp=convWeightImage_p->shape();
675 205 : imshp(2)=npol;
676 :
677 :
678 205 : skyCoverage_p=new TempImage<Float> (imshp, coords,1.0);
679 410 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
680 410 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
681 205 : IPosition stride(4, 1);
682 410 : IPosition trc(blc+image->shape()-stride);
683 :
684 : // Do the copy
685 205 : IPosition start(4, 0);
686 205 : convWeightImage_p->put(griddedWeight(blc, trc));
687 205 : StokesImageUtil::ToStokesPSF(*skyCoverage_p, *convWeightImage_p);
688 205 : 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 205 : pbConvFunc_p->setWeightImage(skyCoverage_p);
706 205 : if(convWeightImage_p) delete convWeightImage_p;
707 205 : convWeightImage_p=0;
708 205 : 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 205 : }
719 :
720 289 : if(!weightLattice.null()) weightLattice=0;
721 289 : griddedWeight.resize();
722 289 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
723 289 : }
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 66601 : void MosaicFT::put(const vi::VisBuffer2& vb, Int row, Bool dopsf,
994 : FTMachine::Type type)
995 : {
996 :
997 :
998 :
999 :
1000 66601 : Timer tim;
1001 66601 : tim.mark();
1002 :
1003 66601 : matchChannel(vb);
1004 :
1005 :
1006 : //cerr << "CHANMAP " << chanMap << endl;
1007 : //No point in reading data if its not matching in frequency
1008 66601 : if(max(chanMap)==-1)
1009 0 : return;
1010 :
1011 : //const Matrix<Float> *imagingweight;
1012 : //imagingweight=&(vb.imagingWeight());
1013 66601 : Matrix<Float> imagingweight;
1014 66601 : getImagingWeight(imagingweight, vb);
1015 :
1016 66601 : if(dopsf) type=FTMachine::PSF;
1017 :
1018 66601 : Cube<Complex> data;
1019 : //Fortran gridder need the flag as ints
1020 66601 : Cube<Int> flags;
1021 66601 : Matrix<Float> elWeight;
1022 66601 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
1023 :
1024 :
1025 :
1026 : Bool iswgtCopy;
1027 : const Float *wgtStorage;
1028 66601 : wgtStorage=elWeight.getStorage(iswgtCopy);
1029 :
1030 :
1031 :
1032 :
1033 : Bool isCopy;
1034 66601 : const Complex *datStorage=0;
1035 :
1036 : // cerr << "dopsf " << dopsf << " isWeightCopy " << iswgtCopy << " " << wgtStorage<< endl;
1037 66601 : if(!dopsf)
1038 40616 : datStorage=data.getStorage(isCopy);
1039 :
1040 :
1041 : // If row is -1 then we pass through all rows
1042 : Int startRow, endRow, nRow;
1043 66601 : if (row==-1) {
1044 66601 : nRow=vb.nRows();
1045 66601 : startRow=0;
1046 66601 : 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 66601 : Matrix<Double> uvw(negateUV(vb));
1056 66601 : Vector<Double> dphase(vb.nRows());
1057 66601 : dphase=0.0;
1058 :
1059 66601 : doUVWRotation_p=true;
1060 66601 : girarUVW(uvw, dphase, vb);
1061 66601 : 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 66601 : 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 66601 : 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 66601 : Int idopsf=0;
1086 66601 : if(dopsf) idopsf=1;
1087 :
1088 :
1089 66601 : Vector<Int> rowFlags(vb.nRows());
1090 66601 : rowFlags=0;
1091 66601 : rowFlags(vb.flagRow())=true;
1092 66601 : if(!usezero_p) {
1093 15772945 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1094 15706344 : 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 66601 : const IPosition& fs=flags.shape();
1126 : //cerr << "flags shape " << fs << endl;
1127 66601 : std::vector<Int>s(fs.begin(), fs.end());
1128 66601 : Int nvp=s[0];
1129 66601 : Int nvc=s[1];
1130 66601 : Int nvisrow=s[2];
1131 66601 : Int csamp=convSampling;
1132 : Bool uvwcopy;
1133 66601 : const Double *uvwstor=uvw.getStorage(uvwcopy);
1134 : Bool gridcopy;
1135 : Bool convcopy;
1136 : Bool wconvcopy;
1137 66601 : const Complex *convstor=convFunc.getStorage(convcopy);
1138 66601 : const Complex *wconvstor=weightConvFunc_p.getStorage(wconvcopy);
1139 66601 : Int nPolConv=convFunc.shape()[2];
1140 66601 : Int nChanConv=convFunc.shape()[3];
1141 66601 : Int nConvFunc=convFunc.shape()(4);
1142 : Bool weightcopy;
1143 : ////////**************************
1144 66601 : Cube<Int> loc(2, nvc, nRow);
1145 66601 : Cube<Int> off(2, nvc, nRow);
1146 66601 : Matrix<Complex> phasor(nvc, nRow);
1147 : Bool delphase;
1148 66601 : Complex * phasorstor=phasor.getStorage(delphase);
1149 66601 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
1150 66601 : const Double * scalestor=uvScale.getStorage(del);
1151 66601 : const Double * offsetstor=uvOffset.getStorage(del);
1152 66601 : Int * locstor=loc.getStorage(del);
1153 66601 : Int * offstor=off.getStorage(del);
1154 66601 : const Double *dpstor=dphase.getStorage(del);
1155 : Int irow;
1156 66601 : Int nth=1;
1157 : #ifdef _OPENMP
1158 66601 : if(numthreads_p >0){
1159 0 : nth=min(numthreads_p, omp_get_max_threads());
1160 : }
1161 : else{
1162 66601 : nth= omp_get_max_threads();
1163 : }
1164 : //nth=min(4,nth);
1165 : #endif
1166 66601 : Double cinv=Double(1.0)/C::c;
1167 :
1168 66601 : Int dow=0;
1169 66601 : #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 66601 : timemass_p +=tim.real();
1184 : Int ixsub, iysub, icounter;
1185 66601 : ixsub=1;
1186 66601 : iysub=1;
1187 : //////***********************DEBUGGING
1188 : //nth=1;
1189 : ////////***************
1190 66601 : if (nth >3){
1191 66601 : ixsub=8;
1192 66601 : iysub=8;
1193 : }
1194 0 : else if(nth >1){
1195 0 : ixsub=2;
1196 0 : iysub=2;
1197 : }
1198 66601 : Int rbeg=startRow+1;
1199 66601 : Int rend=endRow+1;
1200 66601 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
1201 66601 : Vector<Double *> swgtptr(ixsub*iysub);
1202 66601 : Vector<Bool> swgtdel(ixsub*iysub);
1203 4329065 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1204 4262464 : sumwgt[icounter].resize(sumWeight.shape());
1205 4262464 : sumwgt[icounter].set(0.0);
1206 4262464 : swgtptr[icounter]=sumwgt[icounter].getStorage(swgtdel(icounter));
1207 : }
1208 : //cerr << "done thread " << doneThreadPartition_p << " " << ixsub*iysub << endl;
1209 66601 : if(doneThreadPartition_p < 0){
1210 205 : xsect_p.resize(ixsub*iysub);
1211 205 : ysect_p.resize(ixsub*iysub);
1212 205 : nxsect_p.resize(ixsub*iysub);
1213 205 : nysect_p.resize(ixsub*iysub);
1214 13325 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1215 13120 : findGridSector(nx, ny, ixsub, iysub, 0, 0, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
1216 : }
1217 : }
1218 66601 : Vector<Int> xsect, ysect, nxsect, nysect;
1219 66601 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
1220 : //cerr << xsect.shape() << " " << xsect << endl;
1221 66601 : const Int* pmapstor=polMap.getStorage(del);
1222 66601 : const Int* cmapstor=chanMap.getStorage(del);
1223 : // Dummy sumwt for gridweight part
1224 66601 : Matrix<Double> dumSumWeight(npol, nchan);
1225 66601 : dumSumWeight=sumWeight;
1226 : Bool isDSWC;
1227 66601 : Double *dsumwtstor=dumSumWeight.getStorage(isDSWC);
1228 66601 : Int nc=nchan;
1229 66601 : Int np=npol;
1230 66601 : Int nxp=nx;
1231 66601 : Int nyp=ny;
1232 66601 : Int csize=convSize;
1233 66601 : const Int * flagstor=flags.getStorage(del);
1234 66601 : const Int * rowflagstor=rowFlags.getStorage(del);
1235 66601 : Int csupp=convSupport;
1236 66601 : const Int *convrowmapstor=convRowMap_p.getStorage(del);
1237 66601 : const Int *convchanmapstor=convChanMap_p.getStorage(del);
1238 66601 : const Int *convpolmapstor=convPolMap_p.getStorage(del);
1239 : ///
1240 :
1241 :
1242 : ////////***************************
1243 66601 : tim.mark();
1244 :
1245 66601 : if(useDoubleGrid_p) {
1246 66601 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
1247 :
1248 66601 : #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 4329065 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
1302 4262464 : sumwgt[icounter].putStorage(swgtptr[icounter],swgtdel[icounter]);
1303 4262464 : sumWeight=sumWeight+sumwgt[icounter];
1304 : }
1305 :
1306 : //cerr << "SUMWEIG " << sumWeight << endl;
1307 66601 : griddedData2.putStorage(gridstor, gridcopy);
1308 66601 : if(dopsf && (nth >4))
1309 10440 : tweakGridSector(nx, ny, ixsub, iysub);
1310 66601 : timegrid_p+=tim.real();
1311 66601 : tim.mark();
1312 66601 : 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 44379 : DComplex *gridwgtstor=griddedWeight2.getStorage(weightcopy);
1316 44379 : 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 44379 : griddedWeight2.putStorage(gridwgtstor, weightcopy);
1324 :
1325 : }
1326 66601 : 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 66601 : convFunc.freeStorage(convstor, convcopy);
1409 66601 : weightConvFunc_p.freeStorage(wconvstor, wconvcopy);
1410 66601 : dumSumWeight.putStorage(dsumwtstor, isDSWC);
1411 : //cerr << "dumSumwe " << dumSumWeight << endl;
1412 66601 : uvw.freeStorage(uvwstor, uvwcopy);
1413 66601 : if(!dopsf)
1414 40616 : data.freeStorage(datStorage, isCopy);
1415 :
1416 66601 : elWeight.freeStorage(wgtStorage,iswgtCopy);
1417 :
1418 :
1419 :
1420 :
1421 66601 : }
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 86517 : void MosaicFT::ok() {
2253 86517 : AlwaysAssert(image, AipsError);
2254 86517 : }
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 :
|
