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