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