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