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