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