Line data Source code
1 : //# SDGrid.cc: Implementation of SDGrid 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: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Arrays/Array.h>
29 : #include <casacore/casa/Arrays/ArrayLogical.h>
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Arrays/Cube.h>
32 : #include <casacore/casa/Arrays/MaskedArray.h>
33 : #include <casacore/casa/Arrays/Matrix.h>
34 : #include <casacore/casa/Arrays/MatrixIter.h>
35 : #include <casacore/casa/Arrays/Slice.h>
36 : #include <casacore/casa/Arrays/Vector.h>
37 : #include <casacore/casa/BasicSL/Constants.h>
38 : #include <casacore/casa/BasicSL/String.h>
39 : #include <casacore/casa/Containers/Block.h>
40 : #include <casacore/casa/Exceptions/Error.h>
41 : #include <casacore/casa/OS/Timer.h>
42 : #include <casacore/casa/Quanta/MVAngle.h>
43 : #include <casacore/casa/Quanta/MVTime.h>
44 : #include <casacore/casa/Quanta/UnitMap.h>
45 : #include <casacore/casa/Quanta/UnitVal.h>
46 : #include <sstream>
47 : #include <casacore/casa/Utilities/Assert.h>
48 :
49 : #include <components/ComponentModels/ConstantSpectrum.h>
50 : #include <components/ComponentModels/Flux.h>
51 : #include <components/ComponentModels/PointShape.h>
52 :
53 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
54 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
55 : #include <casacore/coordinates/Coordinates/Projection.h>
56 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
57 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
58 :
59 : #include <casacore/images/Images/ImageInterface.h>
60 : #include <casacore/images/Images/PagedImage.h>
61 : #include <casacore/images/Images/TempImage.h>
62 :
63 : #include <casacore/lattices/Lattices/ArrayLattice.h>
64 : #include <casacore/lattices/Lattices/LatticeCache.h>
65 : #include <casacore/lattices/Lattices/LatticeIterator.h>
66 : #include <casacore/lattices/Lattices/LatticeStepper.h>
67 : #include <casacore/lattices/Lattices/SubLattice.h>
68 : #include <casacore/lattices/LEL/LatticeExpr.h>
69 : #include <casacore/lattices/LRegions/LCBox.h>
70 :
71 : #include <casacore/measures/Measures/Stokes.h>
72 : #include <casacore/ms/MeasurementSets/MSColumns.h>
73 : #include <msvis/MSVis/StokesVector.h>
74 : #include <msvis/MSVis/VisBuffer2.h>
75 : #include <msvis/MSVis/VisibilityIterator2.h>
76 : #include <casacore/scimath/Mathematics/RigidVector.h>
77 : #include <synthesis/TransformMachines2/SDGrid.h>
78 : #include <synthesis/TransformMachines2/SkyJones.h>
79 : #include <synthesis/TransformMachines/StokesImageUtil.h>
80 :
81 : #include <casacore/tables/TaQL/TableParse.h>
82 :
83 : using namespace casacore;
84 : namespace casa {
85 : namespace refim {//# namespace for imaging refactor
86 :
87 0 : SDGrid::SDGrid(SkyJones& sj, Int icachesize, Int itilesize,
88 0 : String iconvType, Int userSupport, Bool useImagingWeight)
89 0 : : FTMachine(), sj_p(&sj), imageCache(0), wImageCache(0),
90 0 : cachesize(icachesize), tilesize(itilesize),
91 0 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
92 0 : pointingToImage(0), userSetSupport_p(userSupport),
93 0 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
94 0 : minWeight_p(0.), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
95 0 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(false)
96 : {
97 0 : lastIndex_p=0;
98 0 : }
99 :
100 36 : SDGrid::SDGrid(MPosition& mLocation, SkyJones& sj, Int icachesize, Int itilesize,
101 36 : String iconvType, Int userSupport, Float minweight, Bool clipminmax, Bool useImagingWeight)
102 36 : : FTMachine(), sj_p(&sj), imageCache(0), wImageCache(0),
103 36 : cachesize(icachesize), tilesize(itilesize),
104 36 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
105 72 : pointingToImage(0), userSetSupport_p(userSupport),
106 36 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
107 36 : minWeight_p(minweight), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
108 72 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(clipminmax)
109 : {
110 36 : mLocation_p=mLocation;
111 36 : lastIndex_p=0;
112 36 : }
113 :
114 0 : SDGrid::SDGrid(Int icachesize, Int itilesize,
115 0 : String iconvType, Int userSupport, Bool useImagingWeight)
116 0 : : FTMachine(), sj_p(0), imageCache(0), wImageCache(0),
117 0 : cachesize(icachesize), tilesize(itilesize),
118 0 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
119 0 : pointingToImage(0), userSetSupport_p(userSupport),
120 0 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
121 0 : minWeight_p(0.), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
122 0 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(false)
123 : {
124 0 : lastIndex_p=0;
125 0 : }
126 :
127 232 : SDGrid::SDGrid(MPosition &mLocation, Int icachesize, Int itilesize,
128 232 : String iconvType, Int userSupport, Float minweight, Bool clipminmax, Bool useImagingWeight)
129 232 : : FTMachine(), sj_p(0), imageCache(0), wImageCache(0),
130 232 : cachesize(icachesize), tilesize(itilesize),
131 232 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
132 464 : pointingToImage(0), userSetSupport_p(userSupport),
133 232 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
134 232 : minWeight_p(minweight), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1),
135 232 : msId_p(-1),
136 464 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(clipminmax)
137 : {
138 232 : mLocation_p=mLocation;
139 232 : lastIndex_p=0;
140 232 : }
141 :
142 6 : SDGrid::SDGrid(MPosition &mLocation, Int icachesize, Int itilesize,
143 : String iconvType, Float truncate, Float gwidth, Float jwidth,
144 6 : Float minweight, Bool clipminmax, Bool useImagingWeight)
145 6 : : FTMachine(), sj_p(0), imageCache(0), wImageCache(0),
146 6 : cachesize(icachesize), tilesize(itilesize),
147 6 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
148 12 : pointingToImage(0), userSetSupport_p(-1),
149 6 : truncate_p(truncate), gwidth_p(gwidth), jwidth_p(jwidth),
150 6 : minWeight_p(minweight), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
151 12 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(clipminmax)
152 : {
153 6 : mLocation_p=mLocation;
154 6 : lastIndex_p=0;
155 6 : }
156 :
157 :
158 : //----------------------------------------------------------------------
159 0 : SDGrid& SDGrid::operator=(const SDGrid& other)
160 : {
161 0 : if(this!=&other) {
162 : //Do the base parameters
163 0 : FTMachine::operator=(other);
164 0 : sj_p=other.sj_p;
165 0 : imageCache=other.imageCache;
166 0 : wImage=other.wImage;
167 0 : wImageCache=other.wImageCache;
168 0 : cachesize=other.cachesize;
169 0 : tilesize=other.tilesize;
170 0 : isTiled=other.isTiled;
171 0 : lattice=other.lattice;
172 0 : arrayLattice=other.arrayLattice;
173 0 : wLattice=other.wLattice;
174 0 : wArrayLattice=other.wArrayLattice;
175 0 : convType=other.convType;
176 0 : if(other.wImage !=0)
177 0 : wImage=(TempImage<Float> *)other.wImage->cloneII();
178 : else
179 0 : wImage=0;
180 0 : pointingDirCol_p=other.pointingDirCol_p;
181 0 : pointingToImage=0;
182 0 : xyPos.resize();
183 0 : xyPos=other.xyPos;
184 0 : xyPosMovingOrig_p=other.xyPosMovingOrig_p;
185 0 : convFunc.resize();
186 0 : convFunc=other.convFunc;
187 0 : convSampling=other.convSampling;
188 0 : convSize=other.convSize;
189 0 : convSupport=other.convSupport;
190 0 : userSetSupport_p=other.userSetSupport_p;
191 0 : lastIndex_p=0;
192 0 : lastIndexPerAnt_p=0;
193 0 : lastAntID_p=-1;
194 0 : msId_p=-1;
195 0 : useImagingWeight_p=other.useImagingWeight_p;
196 0 : clipminmax_=other.clipminmax_;
197 : };
198 0 : return *this;
199 : };
200 :
201 548 : String SDGrid::name() const{
202 548 : return String("SDGrid");
203 : }
204 :
205 : //----------------------------------------------------------------------
206 : // Odds are that it changed.....
207 0 : Bool SDGrid::changed(const vi::VisBuffer2& /*vb*/) {
208 0 : return false;
209 : }
210 :
211 : //----------------------------------------------------------------------
212 0 : SDGrid::SDGrid(const SDGrid& other):FTMachine()
213 : {
214 0 : operator=(other);
215 0 : }
216 :
217 : #define NEED_UNDERSCORES
218 : #if defined(NEED_UNDERSCORES)
219 : #define grdsf grdsf_
220 : #define grdgauss grdgauss_
221 : #define grdjinc1 grdjinc1_
222 : #endif
223 :
224 : extern "C" {
225 : void grdsf(Double*, Double*);
226 : void grdgauss(Double*, Double*, Double*);
227 : void grdjinc1(Double*, Double*, Int*, Double*);
228 : }
229 :
230 : //----------------------------------------------------------------------
231 137 : void SDGrid::init() {
232 :
233 : // FIXME: don't mess with parent's class logger
234 : // unless you make sure you reset it to it's original state
235 : // when you are done
236 137 : logIO() << LogOrigin("SDGrid", "init") << LogIO::NORMAL;
237 :
238 : //pfile = fopen("ptdata.txt","w");
239 :
240 137 : ok();
241 :
242 : { // Initialize members
243 137 : isTiled = false;
244 :
245 137 : nx = image->shape()(0);
246 137 : ny = image->shape()(1);
247 137 : npol = image->shape()(2);
248 137 : nchan = image->shape()(3);
249 :
250 137 : sumWeight.resize(npol, nchan);
251 :
252 : // Set up image cache needed for gridding.
253 137 : if (imageCache) delete imageCache;
254 137 : imageCache = 0;
255 :
256 : // Initialize weight image
257 137 : if (wImage) delete wImage;
258 137 : wImage = 0;
259 274 : wImage = new TempImage<Float>(
260 274 : image->shape(),
261 137 : image->coordinates()
262 137 : );
263 : }
264 :
265 : { // Compute Convolution Function
266 137 : convType = downcase(convType);
267 : logIO() << LogIO::NORMAL2
268 137 : << "Convolution function: " << convType
269 137 : << LogIO::POST;
270 :
271 137 : if (convType == "pb") { // Primary Beam: do nothing
272 : //cerr << "CNVFunc " << convFunc << endl;
273 : }
274 119 : else if (convType == "box") { // Box Function
275 99 : convSupport = (userSetSupport_p >= 0) ? userSetSupport_p : 0;
276 : logIO() << LogIO::NORMAL2
277 : << "Support: " << convSupport << " pixels"
278 99 : << LogIO::POST;
279 :
280 99 : convSampling = 100;
281 99 : convSize = convSampling*(2*convSupport+2);
282 99 : convFunc.resize(convSize);
283 99 : convFunc = 0.0;
284 9999 : for (Int i=0; i<convSize/2; i++) {
285 9900 : convFunc(i) = 1.0;
286 : }
287 : }
288 20 : else if (convType == "sf") { // Prolate Spheroidal Wave Function
289 17 : convSupport = (userSetSupport_p >= 0) ? userSetSupport_p : 3;
290 : logIO() << LogIO::NORMAL2
291 : << "Support: " << convSupport << " pixels"
292 17 : << LogIO::POST;
293 :
294 : // FIXME: why 100 ?
295 17 : convSampling = 100;
296 17 : convSize = convSampling*(2*convSupport + 2);
297 17 : convFunc.resize(convSize);
298 17 : convFunc = 0.0;
299 9917 : for (Int i=0; i<convSampling*convSupport; i++) {
300 9900 : Double nu = Double(i)/Double(convSupport*convSampling);
301 : Double val;
302 9900 : grdsf(&nu, &val);
303 9900 : convFunc(i) = (1.0-nu*nu)*val;
304 : }
305 : }
306 3 : else if (convType == "gauss") { // Gauss function
307 : // default is b=1.0 (Mangum et al. 2007)
308 : // FIXME: how does b=1.0 relate to current code ?
309 1 : Double hwhm = (gwidth_p > 0.0) ? Double(gwidth_p) : sqrt(log(2.0));
310 :
311 1 : Float truncate = (truncate_p >= 0.0) ? truncate_p : 3.0 * hwhm;
312 1 : convSampling = 100;
313 1 : Int itruncate = (Int)(truncate*Double(convSampling) + 0.5);
314 :
315 : logIO() << LogIO::NORMAL2
316 : << "hwhm=" << hwhm
317 1 : << LogIO::POST;
318 :
319 : logIO() << LogIO::NORMAL2
320 : << "itruncate=" << itruncate
321 1 : << LogIO::POST;
322 :
323 1 : convSupport = (Int)(truncate);
324 1 : convSupport += ((truncate - (Float)convSupport) > 0.0) ? 1 : 0;
325 :
326 1 : convSize = convSampling*(2*convSupport + 2);
327 :
328 1 : convFunc.resize(convSize);
329 1 : convFunc = 0.0;
330 : Double val, x;
331 252 : for (Int i=0; i<=itruncate; i++) {
332 251 : x = Double(i)/Double(convSampling);
333 251 : grdgauss(&hwhm, &x, &val);
334 251 : convFunc(i)=val;
335 : }
336 :
337 : // String outfile = convType + ".dat";
338 : // ofstream ofs(outfile.c_str());
339 : // for (Int i = 0 ; i < convSize ; i++) {
340 : // ofs << i << " " << convFunc[i] << endl;
341 : // }
342 : // ofs.close();
343 :
344 : }
345 2 : else if (convType == "gjinc") { // Gauss * Jinc function
346 : // default is b=2.52, c=1.55 (Mangum et al. 2007)
347 2 : Double hwhm = (gwidth_p > 0.0) ? Double(gwidth_p) : sqrt(log(2.0))*2.52;
348 2 : Double c = (jwidth_p > 0.0) ? Double(jwidth_p) : 1.55;
349 2 : convSampling = 100;
350 2 : Int itruncate=(Int)(truncate_p*Double(convSampling) + 0.5);
351 :
352 : logIO() << LogIO::NORMAL2
353 : << "hwhm=" << hwhm
354 2 : << LogIO::POST;
355 : logIO() << LogIO::NORMAL2
356 : << "c=" << c
357 2 : << LogIO::POST;
358 : logIO() << LogIO::NORMAL2
359 : << "itruncate=" << itruncate
360 2 : << LogIO::POST;
361 :
362 2 : Float convSupportF = (truncate_p >= 0.0) ? truncate_p : (2*c);
363 2 : convSupport = (Int)convSupportF;
364 2 : convSupport += (((convSupportF-(Float)convSupport) > 0.0) ? 1 : 0);
365 2 : convSize = convSampling*(2*convSupport + 2);
366 :
367 2 : convFunc.resize(convSize);
368 2 : convFunc = 0.0;
369 : Double x, val1, val2;
370 2 : Int normalize = 1;
371 :
372 2 : if (itruncate >= 0) {
373 0 : for (Int i=0 ; i<itruncate; i++) {
374 0 : x = Double(i) / Double(convSampling);
375 0 : grdgauss(&hwhm, &x, &val1);
376 0 : grdjinc1(&c, &x, &normalize, &val2);
377 0 : convFunc(i) = val1 * val2;
378 : }
379 : }
380 : else { // default is to truncate at first null
381 382 : for (Int i=0; i<convSize; i++) {
382 382 : x = Double(i) / Double(convSampling);
383 382 : grdjinc1(&c, &x, &normalize, &val2);
384 382 : if (val2 <= 0.0) {
385 : logIO() << LogIO::NORMAL3
386 : << "convFunc is automatically truncated at radius " << x
387 2 : << LogIO::POST;
388 2 : break;
389 : }
390 380 : grdgauss(&hwhm, &x, &val1);
391 380 : convFunc(i) = val1 * val2;
392 : }
393 : }
394 :
395 : // String outfile = convType + ".dat";
396 : // ofstream ofs(outfile.c_str());
397 : // for (Int i = 0 ; i < convSize ; i++) {
398 : // ofs << i << " " << convFunc[i] << endl;
399 : // }
400 : // ofs.close();
401 :
402 : }
403 : else { // Throw exception
404 0 : logIO_p << "Unknown convolution function: " << convType
405 0 : << LogIO::EXCEPTION;
406 : }
407 : }
408 :
409 137 : }
410 :
411 : // This is nasty, we should use CountedPointers here.
412 548 : SDGrid::~SDGrid() {
413 : //fclose(pfile);
414 274 : if (imageCache) delete imageCache; imageCache = 0;
415 274 : if (arrayLattice) delete arrayLattice; arrayLattice = 0;
416 274 : if (wImage) delete wImage; wImage = 0;
417 274 : if (wImageCache) delete wImageCache; wImageCache = 0;
418 274 : if (wArrayLattice) delete wArrayLattice; wArrayLattice = 0;
419 274 : if (interpolator) delete interpolator; interpolator = 0;
420 548 : }
421 :
422 137 : void SDGrid::initUVWMachine(const vi::VisBuffer2& vb) {
423 : // UVWMachine is not necessary for single dish imaging
424 137 : if (uvwMachine_p) delete uvwMachine_p; uvwMachine_p = 0;
425 137 : phaseShifter_p.reset();
426 :
427 137 : doUVWRotation_p=false;
428 137 : }
429 :
430 18 : void SDGrid::findPBAsConvFunction(const ImageInterface<Complex>& image,
431 : const vi::VisBuffer2& vb) {
432 :
433 : // Get the coordinate system and increase the sampling by
434 : // a factor of ~ 100.
435 18 : CoordinateSystem coords(image.coordinates());
436 :
437 : // Set up the convolution function: make the buffer plenty
438 : // big so that we can trim it back
439 18 : convSupport=max(128, sj_p->support(vb, coords));
440 :
441 18 : convSampling=100;
442 18 : convSize=convSampling*convSupport;
443 :
444 : // Make a one dimensional image to calculate the
445 : // primary beam. We oversample this by a factor of
446 : // convSampling.
447 18 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
448 18 : AlwaysAssert(directionIndex>=0, AipsError);
449 18 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
450 18 : Vector<Double> sampling;
451 18 : sampling = dc.increment();
452 18 : sampling*=1.0/Double(convSampling);
453 18 : dc.setIncrement(sampling);
454 :
455 : // Set the reference value to the first pointing in the coordinate
456 : // system used in the POINTING table.
457 : {
458 18 : uInt row = 0;
459 :
460 : // reset lastAntID_p to use correct antenna position
461 18 : lastAntID_p = -1;
462 :
463 18 : const MSPointingColumns& act_mspc = vb.subtableColumns().pointing();
464 18 : Bool nullPointingTable = (act_mspc.nrow() < 1);
465 18 : Int pointIndex = -1;
466 18 : if (!nullPointingTable){
467 : //if(vb.newMS()) This thing is buggy...using msId change
468 18 : if (vb.msId() != msId_p) {
469 18 : lastIndex_p=0;
470 18 : if (lastIndexPerAnt_p.nelements() < (size_t)vb.nAntennas()) {
471 18 : lastIndexPerAnt_p.resize(vb.nAntennas());
472 : }
473 18 : lastIndexPerAnt_p = 0;
474 18 : msId_p = vb.msId();
475 : }
476 18 : pointIndex = getIndex(act_mspc, vb.time()(row), -1.0, vb.antenna1()(row));
477 18 : if (pointIndex < 0)
478 0 : pointIndex = getIndex(act_mspc, vb.time()(row), vb.timeInterval()(row), vb.antenna1()(row));
479 : }
480 18 : if (!nullPointingTable && ((pointIndex < 0) || (pointIndex >= Int(act_mspc.time().nrow())))) {
481 0 : ostringstream o;
482 0 : o << "Failed to find pointing information for time " <<
483 0 : MVTime(vb.time()(row)/86400.0);
484 0 : logIO_p << String(o) << LogIO::EXCEPTION;
485 0 : }
486 :
487 36 : MEpoch epoch(Quantity(vb.time()(row), "s"));
488 18 : if (!pointingToImage) {
489 18 : lastAntID_p = vb.antenna1()(row);
490 18 : MPosition pos = vb.subtableColumns().antenna().positionMeas()(lastAntID_p);
491 : //mFrame_p = MeasFrame(epoch, pos);
492 18 : (!mFrame_p.epoch()) ? mFrame_p.set(epoch) : mFrame_p.resetEpoch(epoch);
493 18 : (!mFrame_p.position()) ? mFrame_p.set(pos) : mFrame_p.resetPosition(pos);
494 18 : if (!nullPointingTable) {
495 18 : worldPosMeas = directionMeas(act_mspc, pointIndex);
496 : } else {
497 0 : worldPosMeas = vb.direction1()(row);
498 : }
499 :
500 : // Make a machine to convert from the worldPosMeas to the output
501 : // Direction Measure type for the relevant frame
502 18 : MDirection::Ref outRef(dc.directionType(), mFrame_p);
503 18 : pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
504 18 : if (!pointingToImage) {
505 0 : logIO_p << "Cannot make direction conversion machine" << LogIO::EXCEPTION;
506 : }
507 :
508 18 : } else {
509 0 : mFrame_p.resetEpoch(epoch);
510 0 : if (lastAntID_p != vb.antenna1()(row)) {
511 0 : logIO_p << LogIO::DEBUGGING
512 : << "updating antenna position. MS ID " << msId_p
513 : << ", last antenna ID " << lastAntID_p
514 0 : << " new antenna ID " << vb.antenna1()(row) << LogIO::POST;
515 0 : lastAntID_p = vb.antenna1()(row);
516 0 : MPosition pos = vb.subtableColumns().antenna().positionMeas()(lastAntID_p);
517 0 : mFrame_p.resetPosition(pos);
518 0 : }
519 : }
520 :
521 18 : if (!nullPointingTable) {
522 18 : worldPosMeas = (*pointingToImage)(directionMeas(act_mspc, pointIndex));
523 : } else {
524 0 : worldPosMeas = (*pointingToImage)(vb.direction1()(row));
525 : }
526 18 : delete pointingToImage;
527 18 : pointingToImage = 0;
528 18 : }
529 :
530 18 : directionCoord = coords.directionCoordinate(directionIndex);
531 : //make sure we use the same units
532 18 : worldPosMeas.set(dc.worldAxisUnits()(0));
533 :
534 : // Reference pixel may be modified in dc.setReferenceValue when
535 : // projection type is SFL. To take into account this effect,
536 : // keep original reference pixel here and subtract it from
537 : // the reference pixel after dc.setReferenceValue instead
538 : // of setting reference pixel to (0,0).
539 18 : Vector<Double> const originalReferencePixel = dc.referencePixel();
540 18 : dc.setReferenceValue(worldPosMeas.getAngle().getValue());
541 : //Vector<Double> unitVec(2);
542 : //unitVec=0.0;
543 : //dc.setReferencePixel(unitVec);
544 36 : Vector<Double> updatedReferencePixel = dc.referencePixel() - originalReferencePixel;
545 18 : dc.setReferencePixel(updatedReferencePixel);
546 :
547 18 : coords.replaceCoordinate(dc, directionIndex);
548 :
549 18 : IPosition pbShape(4, convSize, 2, 1, 1);
550 18 : IPosition start(4, 0, 0, 0, 0);
551 :
552 18 : TempImage<Complex> onedPB(pbShape, coords);
553 :
554 18 : onedPB.set(Complex(1.0, 0.0));
555 :
556 18 : AlwaysAssert(sj_p, AipsError);
557 18 : sj_p->apply(onedPB, onedPB, vb, 0);
558 :
559 18 : IPosition pbSlice(4, convSize, 1, 1, 1);
560 36 : Vector<Float> tempConvFunc=real(onedPB.getSlice(start, pbSlice, true));
561 : // Find number of significant points
562 18 : uInt cfLen=0;
563 7475 : for(uInt i=0;i<tempConvFunc.nelements();++i) {
564 7475 : if(tempConvFunc(i)<1e-3) break;
565 7457 : ++cfLen;
566 : }
567 18 : if(cfLen<1) {
568 : logIO() << LogIO::SEVERE
569 : << "Possible problem in primary beam calculation: no points in gridding function"
570 0 : << " - no points to be gridded on this image?" << LogIO::POST;
571 0 : cfLen=1;
572 : }
573 18 : Vector<Float> trimConvFunc=tempConvFunc(Slice(0,cfLen-1,1));
574 :
575 : // Now fill in the convolution function vector
576 18 : convSupport=cfLen/convSampling;
577 18 : convSize=convSampling*(2*convSupport+2);
578 18 : convFunc.resize(convSize);
579 18 : convFunc=0.0;
580 18 : convFunc(Slice(0,cfLen-1,1))=trimConvFunc(Slice(0,cfLen-1,1));
581 :
582 :
583 18 : }
584 :
585 : // Initialize for a transform from the Sky domain. This means that
586 : // we grid-correct, and FFT the image
587 0 : void SDGrid::initializeToVis(ImageInterface<Complex>& iimage,
588 : const vi::VisBuffer2& vb)
589 : {
590 0 : image=&iimage;
591 :
592 0 : ok();
593 :
594 0 : init();
595 :
596 0 : if(convType=="pb") {
597 0 : findPBAsConvFunction(*image, vb);
598 : }
599 :
600 : // reset msId_p and lastAntID_p to -1
601 : // this is to ensure correct antenna position in getXYPos
602 0 : msId_p = -1;
603 0 : lastAntID_p = -1;
604 :
605 : // Initialize the maps for polarization and channel. These maps
606 : // translate visibility indices into image indices
607 0 : initMaps(vb);
608 :
609 : // First get the CoordinateSystem for the image and then find
610 : // the DirectionCoordinate
611 0 : CoordinateSystem coords=image->coordinates();
612 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
613 0 : AlwaysAssert(directionIndex>=0, AipsError);
614 0 : directionCoord=coords.directionCoordinate(directionIndex);
615 : /*if((image->shape().product())>cachesize) {
616 : isTiled=true;
617 : }
618 : else {
619 : isTiled=false;
620 : }*/
621 0 : isTiled=false;
622 0 : nx = image->shape()(0);
623 0 : ny = image->shape()(1);
624 0 : npol = image->shape()(2);
625 0 : nchan = image->shape()(3);
626 :
627 : // If we are memory-based then read the image in and create an
628 : // ArrayLattice otherwise just use the PagedImage
629 : /*if(isTiled) {
630 : lattice=image;
631 : wLattice=wImage;
632 : }
633 : else*/
634 : {
635 : // Make the grid the correct shape and turn it into an array lattice
636 0 : IPosition gridShape(4, nx, ny, npol, nchan);
637 0 : griddedData.resize(gridShape);
638 0 : griddedData = Complex(0.0);
639 :
640 0 : wGriddedData.resize(gridShape);
641 0 : wGriddedData = 0.0;
642 :
643 0 : if(arrayLattice) delete arrayLattice; arrayLattice=0;
644 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
645 :
646 0 : if(wArrayLattice) delete wArrayLattice; wArrayLattice=0;
647 0 : wArrayLattice = new ArrayLattice<Float>(wGriddedData);
648 0 : wArrayLattice->set(0.0);
649 0 : wLattice=wArrayLattice;
650 :
651 : // Now find the SubLattice corresponding to the image
652 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
653 0 : IPosition stride(4, 1);
654 0 : IPosition trc(blc+image->shape()-stride);
655 0 : LCBox gridBox(blc, trc, gridShape);
656 0 : SubLattice<Complex> gridSub(*arrayLattice, gridBox, true);
657 :
658 : // Do the copy
659 0 : gridSub.copyData(*image);
660 :
661 0 : lattice=arrayLattice;
662 0 : }
663 0 : AlwaysAssert(lattice, AipsError);
664 0 : AlwaysAssert(wLattice, AipsError);
665 0 : }
666 :
667 0 : void SDGrid::finalizeToVis()
668 : {
669 : /*if(isTiled) {
670 :
671 : logIO() << LogOrigin("SDGrid", "finalizeToVis") << LogIO::NORMAL;
672 :
673 : AlwaysAssert(imageCache, AipsError);
674 : AlwaysAssert(image, AipsError);
675 : ostringstream o;
676 : imageCache->flush();
677 : imageCache->showCacheStatistics(o);
678 : logIO() << o.str() << LogIO::POST;
679 : }*/
680 0 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
681 0 : }
682 :
683 :
684 : // Initialize the FFT to the Sky.
685 : // Here we have to setup and initialize the grid.
686 137 : void SDGrid::initializeToSky(ImageInterface<Complex>& iimage,
687 : Matrix<Float>& weight, const vi::VisBuffer2& vb)
688 : {
689 : // image always points to the image
690 137 : image=&iimage;
691 :
692 137 : ok();
693 :
694 137 : init();
695 :
696 137 : if(convType=="pb") {
697 18 : findPBAsConvFunction(*image, vb);
698 : }
699 :
700 : // reset msId_p and lastAntID_p to -1
701 : // this is to ensure correct antenna position in getXYPos
702 137 : msId_p = -1;
703 137 : lastAntID_p = -1;
704 :
705 : // Initialize the maps for polarization and channel. These maps
706 : // translate visibility indices into image indices
707 137 : initMaps(vb);
708 : //cerr << "ToSky cachesize " << cachesize << " im shape " << (image->shape().product()) << endl;
709 : /*if((image->shape().product())>cachesize) {
710 : isTiled=true;
711 : }
712 : else {
713 : isTiled=false;
714 : }
715 : */
716 : //////////////No longer using isTiled
717 137 : isTiled=false;
718 137 : nx = image->shape()(0);
719 137 : ny = image->shape()(1);
720 137 : npol = image->shape()(2);
721 137 : nchan = image->shape()(3);
722 :
723 137 : sumWeight=0.0;
724 137 : weight.resize(sumWeight.shape());
725 137 : weight=0.0;
726 :
727 : // First get the CoordinateSystem for the image and then find
728 : // the DirectionCoordinate
729 137 : CoordinateSystem coords=image->coordinates();
730 137 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
731 137 : AlwaysAssert(directionIndex>=0, AipsError);
732 137 : directionCoord=coords.directionCoordinate(directionIndex);
733 :
734 : // Initialize for in memory or to disk gridding. lattice will
735 : // point to the appropriate Lattice, either the ArrayLattice for
736 : // in memory gridding or to the image for to disk gridding.
737 : /*if(isTiled) {
738 : imageCache->flush();
739 : image->set(Complex(0.0));
740 : lattice=image;
741 : wLattice=wImage;
742 : }
743 : else*/
744 : {
745 137 : IPosition gridShape(4, nx, ny, npol, nchan);
746 274 : logIO() << LogOrigin("SDGrid", "initializeToSky", WHERE) << LogIO::DEBUGGING
747 274 : << "gridShape = " << gridShape << LogIO::POST;
748 137 : griddedData.resize(gridShape);
749 137 : griddedData=Complex(0.0);
750 137 : if(arrayLattice) delete arrayLattice; arrayLattice=0;
751 137 : arrayLattice = new ArrayLattice<Complex>(griddedData);
752 137 : lattice=arrayLattice;
753 137 : wGriddedData.resize(gridShape);
754 137 : wGriddedData=0.0;
755 137 : if(wArrayLattice) delete wArrayLattice; wArrayLattice=0;
756 137 : wArrayLattice = new ArrayLattice<Float>(wGriddedData);
757 137 : wLattice=wArrayLattice;
758 :
759 : // clipping related stuff
760 137 : if (clipminmax_) {
761 8 : gmin_.resize(gridShape);
762 8 : gmin_ = Complex(FLT_MAX);
763 8 : gmax_.resize(gridShape);
764 8 : gmax_ = Complex(-FLT_MAX);
765 8 : wmin_.resize(gridShape);
766 8 : wmin_ = 0.0f;
767 8 : wmax_.resize(gridShape);
768 8 : wmax_ = 0.0f;
769 8 : npoints_.resize(gridShape);
770 8 : npoints_ = 0;
771 : }
772 137 : }
773 137 : AlwaysAssert(lattice, AipsError);
774 137 : AlwaysAssert(wLattice, AipsError);
775 137 : }
776 :
777 137 : void SDGrid::finalizeToSky()
778 : {
779 :
780 : // Now we flush the cache and report statistics
781 : // For memory based, we don't write anything out yet.
782 : /*if(isTiled) {
783 : logIO() << LogOrigin("SDGrid", "finalizeToSky") << LogIO::NORMAL;
784 :
785 : AlwaysAssert(image, AipsError);
786 : AlwaysAssert(imageCache, AipsError);
787 : imageCache->flush();
788 : ostringstream o;
789 : imageCache->showCacheStatistics(o);
790 : logIO() << o.str() << LogIO::POST;
791 : }
792 : */
793 :
794 137 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
795 137 : }
796 :
797 0 : Array<Complex>* SDGrid::getDataPointer(const IPosition& centerLoc2D,
798 : Bool readonly) {
799 : Array<Complex>* result;
800 : // Is tiled: get tiles and set up offsets
801 0 : centerLoc(0)=centerLoc2D(0);
802 0 : centerLoc(1)=centerLoc2D(1);
803 0 : result=&imageCache->tile(offsetLoc, centerLoc, readonly);
804 0 : return result;
805 : }
806 0 : Array<Float>* SDGrid::getWDataPointer(const IPosition& centerLoc2D,
807 : Bool readonly) {
808 : Array<Float>* result;
809 : // Is tiled: get tiles and set up offsets
810 0 : centerLoc(0)=centerLoc2D(0);
811 0 : centerLoc(1)=centerLoc2D(1);
812 0 : result=&wImageCache->tile(offsetLoc, centerLoc, readonly);
813 0 : return result;
814 : }
815 :
816 : #define NEED_UNDERSCORES
817 : #if defined(NEED_UNDERSCORES)
818 : #define ggridsd ggridsd_
819 : #define dgridsd dgridsd_
820 : #define ggridsdclip ggridsdclip_
821 : #endif
822 :
823 : extern "C" { // Gridders interfaces
824 : void ggridsd(
825 : Double*,
826 : const Complex*,
827 : Int*,
828 : Int*,
829 : Int*,
830 : const Int*,
831 : const Int*,
832 : const Float*,
833 : Int*,
834 : Int*,
835 : Complex*,
836 : Float*,
837 : Int*,
838 : Int*,
839 : Int*,
840 : Int*,
841 : Int*,
842 : Int*,
843 : Float*,
844 : Int*,
845 : Int*,
846 : Double*
847 : );
848 : void ggridsdclip(
849 : Double*,
850 : const Complex*,
851 : Int*,
852 : Int*,
853 : Int*,
854 : const Int*,
855 : const Int*,
856 : const Float*,
857 : Int*,
858 : Int*,
859 : Complex*,
860 : Float*,
861 : Int*,
862 : Complex*,
863 : Float*,
864 : Complex*,
865 : Float*,
866 : Int*,
867 : Int*,
868 : Int *,
869 : Int *,
870 : Int*,
871 : Int*,
872 : Float*,
873 : Int*,
874 : Int*,
875 : Double*
876 : );
877 : void dgridsd(
878 : Double*,
879 : Complex*,
880 : Int*,
881 : Int*,
882 : const Int*,
883 : const Int*,
884 : Int*,
885 : Int*,
886 : const Complex*,
887 : Int*,
888 : Int*,
889 : Int *,
890 : Int *,
891 : Int*,
892 : Int*,
893 : Float*,
894 : Int*,
895 : Int*
896 : );
897 : }
898 :
899 188441 : void SDGrid::put(const vi::VisBuffer2& vb, Int row, Bool dopsf,
900 : FTMachine::Type type)
901 : {
902 376882 : LogIO os(LogOrigin("SDGrid", "put"));
903 :
904 188441 : gridOk(convSupport);
905 :
906 : // There is no channel mapping cache in VI/VB2 version of FTMachine
907 : // Perform matchChannel everytime
908 188441 : matchChannel(vb);
909 :
910 : // No point in reading data if its not matching in frequency
911 188441 : if (max(chanMap)==-1) return;
912 :
913 169226 : Matrix<Float> imagingweight;
914 : //imagingweight=&(vb.imagingWeight());
915 169226 : pickWeights(vb, imagingweight);
916 :
917 169226 : if (type==FTMachine::PSF || type==FTMachine::COVERAGE) dopsf=true;
918 169226 : if (dopsf) type=FTMachine::PSF;
919 169226 : Cube<Complex> data;
920 169226 : Cube<Int> flags; //Fortran gridder need the flag as ints
921 169226 : Matrix<Float> elWeight;
922 169226 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
923 : //cerr << "number of rows " << vb.nRow()
924 : // << " data shape " << data.shape() << endl;
925 : Bool iswgtCopy;
926 : const Float *wgtStorage;
927 169226 : wgtStorage=elWeight.getStorage(iswgtCopy);
928 : Bool isCopy;
929 169226 : const Complex *datStorage=0;
930 169226 : if (!dopsf) datStorage=data.getStorage(isCopy);
931 :
932 : // If row is -1 then we pass through all rows
933 : Int startRow, endRow, nRow;
934 169226 : if (row==-1) {
935 169226 : nRow=vb.nRows();
936 169226 : startRow=0;
937 169226 : endRow=nRow-1;
938 : } else {
939 0 : nRow=1;
940 0 : startRow=row;
941 0 : endRow=row;
942 : }
943 :
944 169226 : Vector<Int> rowFlags(vb.flagRow().nelements(), 0);
945 345737 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
946 176511 : if(vb.flagRow()(rownr)) rowFlags(rownr)=1;
947 : }
948 :
949 : // Take care of translation of Bools to Integer
950 169226 : Int idopsf = dopsf ? 1 : 0;
951 :
952 : { // Compute spectra pixel coordinates and call gridder
953 : // Make sure failed getXYPos does not fall on grid
954 169226 : constexpr Double kFarAway = -1e9;
955 169226 : Matrix<Double> xyPositions(2, endRow-startRow+1, kFarAway);
956 345737 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
957 176511 : if (getXYPos(vb, rownr)) {
958 176511 : xyPositions(0, rownr)=xyPos(0);
959 176511 : xyPositions(1, rownr)=xyPos(1);
960 : }
961 : }
962 : { // Call gridder
963 : Bool del;
964 169226 : const IPosition& fs=flags.shape();
965 169226 : std::vector<Int> s(fs.begin(), fs.end());
966 : Bool datCopy, wgtCopy;
967 169226 : Complex * datStor=griddedData.getStorage(datCopy);
968 169226 : Float * wgtStor=wGriddedData.getStorage(wgtCopy);
969 :
970 : //Bool call_ggridsd = !clipminmax_ || dopsf;
971 169226 : Bool call_ggridsd = !clipminmax_;
972 :
973 169226 : if (call_ggridsd) { // Call plain gridder
974 :
975 338402 : ggridsd(
976 : xyPositions.getStorage(del),
977 : datStorage,
978 169201 : &s[0],
979 169201 : &s[1],
980 : &idopsf,
981 169201 : flags.getStorage(del),
982 169201 : rowFlags.getStorage(del),
983 : wgtStorage,
984 169201 : &s[2],
985 : &row,
986 : datStor,
987 : wgtStor,
988 : &nx,
989 : &ny,
990 : &npol,
991 : &nchan,
992 : &convSupport,
993 : &convSampling,
994 : convFunc.getStorage(del),
995 : chanMap.getStorage(del),
996 : polMap.getStorage(del),
997 : sumWeight.getStorage(del)
998 : );
999 :
1000 : } else { // Call clipping gridder
1001 : Bool gminCopy;
1002 25 : Complex *gminStor = gmin_.getStorage(gminCopy);
1003 : Bool gmaxCopy;
1004 25 : Complex *gmaxStor = gmax_.getStorage(gmaxCopy);
1005 : Bool wminCopy;
1006 25 : Float *wminStor = wmin_.getStorage(wminCopy);
1007 : Bool wmaxCopy;
1008 25 : Float *wmaxStor = wmax_.getStorage(wmaxCopy);
1009 : Bool npCopy;
1010 25 : Int *npStor = npoints_.getStorage(npCopy);
1011 :
1012 50 : ggridsdclip(
1013 : xyPositions.getStorage(del),
1014 : datStorage,
1015 25 : &s[0],
1016 25 : &s[1],
1017 : &idopsf,
1018 25 : flags.getStorage(del),
1019 25 : rowFlags.getStorage(del),
1020 : wgtStorage,
1021 25 : &s[2],
1022 : &row,
1023 : datStor,
1024 : wgtStor,
1025 : npStor,
1026 : gminStor,
1027 : wminStor,
1028 : gmaxStor,
1029 : wmaxStor,
1030 : &nx,
1031 : &ny,
1032 : &npol,
1033 : &nchan,
1034 : &convSupport,
1035 : &convSampling,
1036 : convFunc.getStorage(del),
1037 : chanMap.getStorage(del),
1038 : polMap.getStorage(del),
1039 : sumWeight.getStorage(del)
1040 : );
1041 :
1042 25 : gmin_.putStorage(gminStor, gminCopy);
1043 25 : gmax_.putStorage(gmaxStor, gmaxCopy);
1044 25 : wmin_.putStorage(wminStor, wminCopy);
1045 25 : wmax_.putStorage(wmaxStor, wmaxCopy);
1046 25 : npoints_.putStorage(npStor, npCopy);
1047 : }
1048 169226 : griddedData.putStorage(datStor, datCopy);
1049 169226 : wGriddedData.putStorage(wgtStor, wgtCopy);
1050 169226 : }
1051 169226 : }
1052 :
1053 : { // Free memory
1054 169226 : if (!dopsf) data.freeStorage(datStorage, isCopy);
1055 :
1056 169226 : elWeight.freeStorage(wgtStorage, iswgtCopy);
1057 : }
1058 :
1059 188441 : }
1060 :
1061 0 : void SDGrid::get(vi::VisBuffer2& vb, Int row)
1062 : {
1063 0 : LogIO os(LogOrigin("SDGrid", "get"));
1064 :
1065 0 : gridOk(convSupport);
1066 :
1067 : // If row is -1 then we pass through all rows
1068 : Int startRow, endRow, nRow;
1069 0 : if (row==-1) {
1070 0 : nRow=vb.nRows();
1071 0 : startRow=0;
1072 0 : endRow=nRow-1;
1073 : // TODO: ask imager guru if commenting out the following line
1074 : // is safe for SDGrid
1075 : //unnecessary zeroing
1076 : //vb.modelVisCube()=Complex(0.0,0.0);
1077 : } else {
1078 0 : nRow=1;
1079 0 : startRow=row;
1080 0 : endRow=row;
1081 : // TODO: ask imager guru if commenting out the following line
1082 : // is safe for SDGrid
1083 : //unnecessary zeroing
1084 : //vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
1085 : }
1086 :
1087 : // There is no channel mapping cache in VI/VB2 version of FTMachine
1088 : // Perform matchChannel everytime
1089 0 : matchChannel(vb);
1090 :
1091 : //No point in reading data if its not matching in frequency
1092 0 : if(max(chanMap)==-1)
1093 0 : return;
1094 0 : Cube<Complex> data;
1095 0 : Cube<Int> flags;
1096 0 : getInterpolateArrays(vb, data, flags);
1097 :
1098 : Complex *datStorage;
1099 : Bool isCopy;
1100 0 : datStorage=data.getStorage(isCopy);
1101 : // NOTE: with MS V2.0 the pointing could change per antenna and timeslot
1102 : //
1103 0 : Vector<Int> rowFlags(vb.flagRow().nelements());
1104 0 : rowFlags=0;
1105 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1106 0 : if(vb.flagRow()(rownr)) rowFlags(rownr)=1;
1107 : //single dish yes ?
1108 0 : if(max(vb.uvw().column(rownr)) > 0.0) rowFlags(rownr)=1;
1109 : }
1110 :
1111 :
1112 : /*if(isTiled) {
1113 :
1114 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1115 :
1116 : if(getXYPos(vb, rownr)) {
1117 :
1118 : // Get the tile
1119 : IPosition centerLoc2D(2, Int(xyPos(0)), Int(xyPos(1)));
1120 : Array<Complex>* dataPtr=getDataPointer(centerLoc2D, true);
1121 : Int aNx=dataPtr->shape()(0);
1122 : Int aNy=dataPtr->shape()(1);
1123 :
1124 : // Now use FORTRAN to do the gridding. Remember to
1125 : // ensure that the shape and offsets of the tile are
1126 : // accounted for.
1127 : Bool del;
1128 : Vector<Double> actualPos(2);
1129 : for (Int i=0;i<2;i++) {
1130 : actualPos(i)=xyPos(i)-Double(offsetLoc(i));
1131 : }
1132 : // IPosition s(data.shape());
1133 : const IPosition& fs=data.shape();
1134 : std::vector<Int> s(fs.begin(), fs.end());
1135 :
1136 : dgridsd(actualPos.getStorage(del),
1137 : datStorage,
1138 : &s[0],
1139 : &s[1],
1140 : flags.getStorage(del),
1141 : rowFlags.getStorage(del),
1142 : &s[2],
1143 : &rownr,
1144 : dataPtr->getStorage(del),
1145 : &aNx,
1146 : &aNy,
1147 : &npol,
1148 : &nchan,
1149 : &convSupport,
1150 : &convSampling,
1151 : convFunc.getStorage(del),
1152 : chanMap.getStorage(del),
1153 : polMap.getStorage(del));
1154 : }
1155 : }
1156 : }
1157 : else*/
1158 : {
1159 0 : Matrix<Double> xyPositions(2, endRow-startRow+1);
1160 0 : xyPositions=-1e9;
1161 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1162 0 : if(getXYPos(vb, rownr)) {
1163 0 : xyPositions(0, rownr)=xyPos(0);
1164 0 : xyPositions(1, rownr)=xyPos(1);
1165 : }
1166 : }
1167 :
1168 : Bool del;
1169 : // IPosition s(data.shape());
1170 0 : const IPosition& fs=data.shape();
1171 0 : std::vector<Int> s(fs.begin(), fs.end());
1172 0 : dgridsd(xyPositions.getStorage(del),
1173 : datStorage,
1174 0 : &s[0],
1175 0 : &s[1],
1176 0 : flags.getStorage(del),
1177 0 : rowFlags.getStorage(del),
1178 0 : &s[2],
1179 : &row,
1180 0 : griddedData.getStorage(del),
1181 : &nx,
1182 : &ny,
1183 : &npol,
1184 : &nchan,
1185 : &convSupport,
1186 : &convSampling,
1187 : convFunc.getStorage(del),
1188 : chanMap.getStorage(del),
1189 : polMap.getStorage(del));
1190 :
1191 0 : data.putStorage(datStorage, isCopy);
1192 0 : }
1193 0 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1194 0 : }
1195 :
1196 :
1197 160 : Bool SDGrid::mustConvertPointingColumn(const MeasurementSet &ms)
1198 : {
1199 160 : const auto havePointings = ms.pointing().nrow() > 0;
1200 160 : if (not havePointings) return false;
1201 :
1202 160 : switch(convertFirst) {
1203 160 : case ConvertFirst::ALWAYS: return true;
1204 160 : case ConvertFirst::NEVER: return false;
1205 0 : case ConvertFirst::AUTO:
1206 : {
1207 0 : const auto nPointings = ms.pointing().nrow();
1208 0 : const auto nSelectedDataRows = ms.nrow();
1209 0 : return nSelectedDataRows > nPointings ? true : false;
1210 : }
1211 0 : default:
1212 0 : LogIO logger(LogOrigin("SDGrid", "mustConvertPointingColumn", WHERE));
1213 : logger << "Bug ! Got unexpected value for SDGrid::convertFirst: "
1214 0 : << int(convertFirst)
1215 0 : << LogIO::EXCEPTION;
1216 : }
1217 0 : return false;
1218 : }
1219 :
1220 0 : void SDGrid::convertPointingColumn(
1221 : const MeasurementSet &ms,
1222 : const MSPointingEnums::PredefinedColumns columnToConvert,
1223 : const MDirection::Types directionRef)
1224 : {
1225 0 : LogIO logger(LogOrigin("SDGrid", "convertPointingColumn"));
1226 :
1227 : const auto & nameOfColumnToConvert =
1228 0 : MSPointing::columnName(columnToConvert);
1229 :
1230 0 : const auto & nameOfDirectionRef = MDirection::showType(directionRef);
1231 :
1232 : logger << "Converting POINTING table column: " << nameOfColumnToConvert
1233 : << " to: " << nameOfDirectionRef
1234 0 : << LogIO::POST;
1235 :
1236 : { // Check parameters
1237 : using POINTING = MSPointingEnums::PredefinedColumns;
1238 : // Column must be a direction column
1239 0 : const auto isDirectionColumn = (
1240 : columnToConvert == POINTING::DIRECTION
1241 0 : or columnToConvert == POINTING::TARGET
1242 0 : or columnToConvert == POINTING::POINTING_OFFSET
1243 0 : or columnToConvert == POINTING::SOURCE_OFFSET
1244 0 : or columnToConvert == POINTING::ENCODER
1245 : );
1246 0 : if (not isDirectionColumn) {
1247 : logger << nameOfColumnToConvert << ": not a direction column"
1248 0 : << LogIO::EXCEPTION;
1249 : }
1250 : }
1251 :
1252 : { // Copy Pointing table structure
1253 0 : constexpr auto doNotCopyRows = True;
1254 0 : ramPointingTable = ms.pointing().copyToMemoryTable(
1255 0 : ms.pointing().tableName() +
1256 0 : "." + MSPointing::columnName(columnToConvert) +
1257 0 : "." + MDirection::showType(directionRef),
1258 : doNotCopyRows
1259 0 : );
1260 0 : ramPointingColumnsPtr.reset(new MSPointingColumns {ramPointingTable});
1261 : }
1262 :
1263 : { // Set the reference frame of the direction columns
1264 : // ---- All direction columns, except the encoder column
1265 0 : ramPointingColumnsPtr->setDirectionRef(directionRef);
1266 : // ---- Encoder column
1267 0 : ramPointingColumnsPtr->setEncoderDirectionRef(directionRef);
1268 : }
1269 :
1270 0 : auto quote = [](const String & s) {
1271 0 : return String("\"") + s + String("\"");
1272 : };
1273 :
1274 : { // Perform 1 dummy conversion:
1275 : // Convert the direction of the first pointing
1276 : // 1 day before it was actually recorded
1277 : // to pre-set static variables in casacore functions like dUT1
1278 : // so that we are sure they will be updated when we convert the column
1279 0 : stringstream dummyConversion; // TaQL command
1280 : { // Create it
1281 : dummyConversion <<
1282 : "using style python\n"
1283 : "select\n"
1284 : " [\n"
1285 : " meas.direction(\n" <<
1286 0 : " " << quote(nameOfDirectionRef) << "\n"
1287 : " , pointing.DIRECTION\n"
1288 : " , (pointing.TIME - 1d), 'UTC'\n"
1289 : " , antenna.POSITION\n"
1290 : " )\n"
1291 : " ] as CONVERTED_DIRECTION_ONE_DAY_BEFORE\n"
1292 : " , (pointing.TIME - 1d) d as oneDayBefore\n"
1293 : " , pointing.TIME d as pointingDay\n"
1294 : " , cdatetime(pointing.TIME) as pointingDay_Str\n"
1295 : " , cdatetime(pointing.TIME -1d) as oneDayBefore_Str\n"
1296 : "from\n"
1297 : " $1 as pointing\n"
1298 : "join\n"
1299 : " $2 as antenna\n"
1300 : " on pointing.ANTENNA_ID == antenna.rowid()\n"
1301 : "where\n"
1302 0 : " pointing.rowid() == 0\n"
1303 : ;
1304 : }
1305 : { // Execute it
1306 : vector<const Table*> tables {
1307 0 : &ms.pointing(), // $1
1308 0 : &ms.antenna(), // $2
1309 0 : };
1310 0 : tableCommand(dummyConversion.str(), tables);
1311 0 : }
1312 0 : }
1313 :
1314 : { // Now convert the column
1315 0 : stringstream convertColumn; // TaQL command
1316 : { // Create it
1317 : convertColumn <<
1318 : "using style python\n"
1319 : "insert\n"
1320 : " into $3 as ram_pointing_table\n"
1321 : " (\n"
1322 : " ANTENNA_ID\n"
1323 : " , TIME\n"
1324 : " , INTERVAL\n"
1325 : " , NUM_POLY\n"
1326 0 : " , " << nameOfColumnToConvert << "\n"
1327 : " )\n"
1328 : "select\n"
1329 : " ANTENNA_ID as ANTENNA_ID_INT INTEGER\n"
1330 : " , TIME\n"
1331 : " , INTERVAL\n"
1332 0 : " , 0 as NUM_POLY INTEGER\n";
1333 : using Pointing = MSPointingEnums::PredefinedColumns;
1334 0 : switch(columnToConvert) {
1335 0 : case Pointing::ENCODER: { // ScalarColumn
1336 : convertColumn <<
1337 : " , meas.direction(\n"
1338 0 : " " << quote(nameOfDirectionRef) << "\n"
1339 0 : " , pointing." << nameOfColumnToConvert << "\n"
1340 : " , pointing.TIME\n"
1341 : " , antenna.POSITION\n"
1342 0 : " )\n"
1343 : ;
1344 0 : break;
1345 : }
1346 0 : default: { // All other direction columns are ArrayColumns
1347 : convertColumn <<
1348 : " , [\n"
1349 : " meas.direction(\n"
1350 0 : " " << quote(nameOfDirectionRef) << "\n"
1351 0 : " , pointing." << nameOfColumnToConvert << "\n"
1352 : " , pointing.TIME\n"
1353 : " , antenna.POSITION\n"
1354 : " )\n"
1355 0 : " ]\n"
1356 : ;
1357 : }
1358 : }
1359 : convertColumn <<
1360 : "from\n"
1361 : " $1 as pointing\n"
1362 : "join\n"
1363 : " $2 as antenna\n"
1364 0 : " on pointing.ANTENNA_ID = antenna.rowid()\n";
1365 : }
1366 : { // Execute it
1367 : vector<const casacore::Table*> tables {
1368 0 : &ms.pointing(), // $1
1369 0 : &ms.antenna(), // $2
1370 : &ramPointingTable // $3
1371 0 : };
1372 0 : tableCommand(convertColumn.str(), tables);
1373 0 : }
1374 0 : }
1375 : logger << "Converted POINTING table column: " << nameOfColumnToConvert
1376 : << " to: " << nameOfDirectionRef
1377 0 : << LogIO::POST;
1378 0 : }
1379 :
1380 0 : void SDGrid::handleNewMs(
1381 : const MeasurementSet &ms,
1382 : ImageInterface<Complex>& image)
1383 : {
1384 0 : if (mustConvertPointingColumn(ms)) {
1385 0 : const auto columnEnum = MSPointing::columnType(pointingDirCol_p);
1386 : const auto imageDirectionRef =
1387 0 : image.coordinates().directionCoordinate().directionType();
1388 0 : convertPointingColumn(ms, columnEnum, imageDirectionRef);
1389 : }
1390 : else {
1391 0 : ramPointingTable = MSPointing();
1392 0 : ramPointingColumnsPtr.reset();
1393 : }
1394 0 : }
1395 :
1396 160 : void SDGrid::handleNewMs(const MeasurementSet & ms,
1397 : CountedPtr<SIImageStore> imstore)
1398 : {
1399 160 : if (imstore.null()) return;
1400 :
1401 160 : if (mustConvertPointingColumn(ms)) {
1402 0 : const auto coordinateSystem = imstore->getCSys();
1403 : const auto imagesDirectionRef =
1404 0 : coordinateSystem.directionCoordinate()
1405 0 : .directionType();
1406 0 : const auto columnEnum = MSPointing::columnType(pointingDirCol_p);
1407 0 : convertPointingColumn(ms, columnEnum, imagesDirectionRef);
1408 0 : }
1409 : else {
1410 160 : ramPointingTable = MSPointing();
1411 160 : ramPointingColumnsPtr.reset();
1412 : }
1413 : }
1414 :
1415 : // Make a plain straightforward honest-to-FSM image. This returns
1416 : // a complex image, without conversion to Stokes. The representation
1417 : // is that required for the visibilities.
1418 : //----------------------------------------------------------------------
1419 0 : void SDGrid::makeImage(
1420 : FTMachine::Type type,
1421 : vi::VisibilityIterator2& vi,
1422 : ImageInterface<Complex>& theImage,
1423 : Matrix<Float>& weight) {
1424 :
1425 0 : logIO() << LogOrigin("FTMachine", "makeImage0") << LogIO::NORMAL;
1426 :
1427 0 : vi::VisBuffer2 *vb = vi.getVisBuffer();
1428 0 : vi.origin();
1429 :
1430 : { // Set Stokes Representation
1431 0 : if (vb->polarizationFrame()==MSIter::Linear) {
1432 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
1433 : }
1434 : else {
1435 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
1436 : }
1437 : }
1438 :
1439 0 : if (type==FTMachine::CORRECTED) { // What if we have no correctedData ?
1440 : const auto haveCorrectedData = not (
1441 0 : MSMainColumns(vi.ms()).correctedData().isNull()
1442 0 : );
1443 0 : if (not haveCorrectedData) {
1444 0 : type = FTMachine::OBSERVED;
1445 : }
1446 : }
1447 :
1448 0 : Bool normalize = (type==FTMachine::COVERAGE) ? false : true;
1449 :
1450 0 : Int Nx = theImage.shape()(0);
1451 0 : Int Ny = theImage.shape()(1);
1452 0 : Int Npol = theImage.shape()(2);
1453 0 : Int Nchan = theImage.shape()(3);
1454 :
1455 0 : Double memtot = Double(HostInfo::memoryTotal(true))*1024.0; // return in kB
1456 0 : Int nchanInMem = Int(memtot/2.0/8.0/Double(Nx*Ny*Npol));
1457 0 : Int nloop = nchanInMem >= Nchan ? 1 : Nchan/nchanInMem+1;
1458 :
1459 0 : ImageInterface<Complex> *imCopy = NULL;
1460 : { // Initialize imCopy if needed
1461 0 : if (nloop==1) {
1462 0 : imCopy = &theImage;
1463 0 : nchanInMem = Nchan;
1464 : }
1465 : else {
1466 : logIO() << "Not enough memory to image in one go \n"
1467 : << " will process the image in " << nloop << " sections"
1468 0 : << LogIO::POST;
1469 : }
1470 : }
1471 :
1472 0 : weight.resize(Npol, Nchan);
1473 0 : Matrix<Float> wgtcopy(Npol, Nchan);
1474 :
1475 0 : Bool isWgtZero = true;
1476 0 : IPosition blc(theImage.shape().size(), 0);
1477 0 : IPosition trc(theImage.shape() - 1);
1478 0 : for (Int k=0; k < nloop; ++k) {
1479 : Int bchan; // Slice boundaries along the channel axis
1480 : Int echan;
1481 : { // Compute them
1482 0 : bchan = k*nchanInMem;
1483 0 : echan = (k+1)*nchanInMem < Nchan ? (k+1)*nchanInMem-1 : Nchan-1;
1484 : }
1485 :
1486 0 : if (nloop > 1) { // Slide. Copy of a slice of theImage
1487 0 : blc[3] = bchan;
1488 0 : trc[3] = echan;
1489 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1490 0 : imCopy = new SubImage<Complex>(theImage, sl, true);
1491 0 : wgtcopy.resize(npol, echan-bchan+1);
1492 0 : }
1493 :
1494 : { // Rewind iterator, initializeToSky
1495 0 : vi.originChunks();
1496 0 : vi.origin();
1497 0 : initializeToSky(*imCopy, wgtcopy, *vb);
1498 : }
1499 :
1500 : { // Debug messages for minmax clipping
1501 0 : logIO() << LogOrigin("SDGrid", "makeImage", WHERE) << LogIO::DEBUGGING
1502 0 : << "doclip_ = " << (clipminmax_ ? "TRUE" : "FALSE")
1503 0 : << " (" << clipminmax_ << ")"
1504 0 : << LogIO::POST;
1505 0 : if (clipminmax_) {
1506 0 : logIO() << LogOrigin("SDGrid", "makeImage", WHERE) << LogIO::DEBUGGING
1507 : << "use ggridsd2 for imaging"
1508 0 : << LogIO::POST;
1509 : }
1510 : }
1511 :
1512 : // Loop over the visibilities, putting VisBuffers
1513 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
1514 0 : if (vi.getImpl()->isNewMs()) {
1515 : // When we pre-convert the user-specified POINTING column
1516 : // - e.g. when convertFirst = always -, re-converting it
1517 : // at each slice iteration is an implementation decision.
1518 : // The slice loop is probably rarely used, and when it is
1519 : // it means we have little RAM available.
1520 0 : handleNewMs(vi.ms(), theImage);
1521 : }
1522 0 : for (vi.origin(); vi.more(); vi.next()) {
1523 0 : switch(type) {
1524 0 : case FTMachine::RESIDUAL:
1525 0 : vb->setVisCube(vb->visCubeCorrected() - vb->visCubeModel());
1526 0 : put(*vb, -1, false);
1527 0 : break;
1528 0 : case FTMachine::MODEL:
1529 0 : put(*vb, -1, false, FTMachine::MODEL);
1530 0 : break;
1531 0 : case FTMachine::CORRECTED:
1532 0 : put(*vb, -1, false, FTMachine::CORRECTED);
1533 0 : break;
1534 0 : case FTMachine::PSF:
1535 0 : vb->setVisCube(Complex(1.0,0.0));
1536 0 : put(*vb, -1, true, FTMachine::PSF);
1537 0 : break;
1538 0 : case FTMachine::COVERAGE:
1539 0 : vb->setVisCube(Complex(1.0));
1540 0 : put(*vb, -1, true, FTMachine::COVERAGE);
1541 0 : break;
1542 0 : case FTMachine::OBSERVED:
1543 : default:
1544 0 : put(*vb, -1, false, FTMachine::OBSERVED);
1545 0 : break;
1546 : }
1547 : }
1548 : }
1549 :
1550 0 : finalizeToSky();
1551 :
1552 : // Normalize by dividing out weights, etc.
1553 0 : getImage(wgtcopy, normalize);
1554 :
1555 : { // Check if all weights are zero
1556 0 : if (max(wgtcopy)==0.0) {
1557 0 : if (nloop > 1) {
1558 : logIO() << LogIO::WARN
1559 : << "No useful data in SDGrid: weights all zero for image slice "
1560 : << k
1561 0 : << LogIO::POST;
1562 : }
1563 : }
1564 : else {
1565 0 : isWgtZero = false;
1566 : }
1567 : }
1568 :
1569 0 : weight(Slice(0, Npol), Slice(bchan, echan-bchan+1)) = wgtcopy;
1570 0 : if (nloop > 1) delete imCopy;
1571 :
1572 : } // loop k
1573 :
1574 0 : if (isWgtZero) { // Log severe warning but don't abort
1575 : logIO() << LogIO::SEVERE
1576 : << "No useful data in SDGrid: weights all zero"
1577 0 : << LogIO::POST;
1578 : }
1579 0 : }
1580 :
1581 :
1582 :
1583 : // Finalize : optionally normalize by weight image
1584 137 : ImageInterface<Complex>& SDGrid::getImage(Matrix<Float>& weights,
1585 : Bool normalize)
1586 : {
1587 137 : AlwaysAssert(lattice, AipsError);
1588 137 : AlwaysAssert(image, AipsError);
1589 :
1590 137 : logIO() << LogOrigin("SDGrid", "getImage") << LogIO::NORMAL;
1591 :
1592 : // execute minmax clipping
1593 137 : clipMinMax();
1594 :
1595 137 : weights.resize(sumWeight.shape());
1596 :
1597 137 : convertArray(weights,sumWeight);
1598 :
1599 : // If the weights are all zero then we cannot normalize
1600 : // otherwise we don't care.
1601 : ///////////////////////
1602 : /*{
1603 : PagedImage<Float> thisScreen(lattice->shape(), image->coordinates(), "TheData");
1604 : LatticeExpr<Float> le(abs(*lattice));
1605 : thisScreen.copyData(le);
1606 : PagedImage<Float> thisScreen2(lattice->shape(), image->coordinates(), "TheWeight");
1607 : LatticeExpr<Float> le2(abs(*wLattice));
1608 : thisScreen2.copyData(le2);
1609 : }*/
1610 : /////////////////////
1611 :
1612 137 : if(normalize) {
1613 0 : if(max(weights)==0.0) {
1614 : //logIO() << LogIO::WARN << "No useful data in SDGrid: weights all zero"
1615 : // << LogIO::POST;
1616 : //image->set(Complex(0.0));
1617 0 : return *image;
1618 : }
1619 : //Timer tim;
1620 : //tim.mark();
1621 : //lattice->copyData((LatticeExpr<Complex>)(iif((*wLattice<=minWeight_p), 0.0,
1622 : // (*lattice)/(*wLattice))));
1623 : //As we are not using disk based lattices anymore the above uses too much memory for
1624 : // ArrayLattice...it does not do a real inplace math but rather mutiple copies
1625 : // seem to be involved thus the less elegant but faster and less memory hog loop
1626 : // below.
1627 :
1628 0 : IPosition pos(4);
1629 0 : for (Int chan=0; chan < nchan; ++chan){
1630 0 : pos[3]=chan;
1631 0 : for( Int pol=0; pol < npol; ++ pol){
1632 0 : pos[2]=pol;
1633 0 : for (Int y=0; y < ny ; ++y){
1634 0 : pos[1]=y;
1635 0 : for (Int x=0; x < nx; ++x){
1636 0 : pos[0]=x;
1637 0 : Float wgt=wGriddedData(pos);
1638 0 : if(wgt > minWeight_p)
1639 0 : griddedData(pos)=griddedData(pos)/wgt;
1640 : else
1641 0 : griddedData(pos)=0.0;
1642 : }
1643 : }
1644 : }
1645 : }
1646 : //tim.show("After normalizing");
1647 0 : }
1648 :
1649 : //if(!isTiled)
1650 : {
1651 : // Now find the SubLattice corresponding to the image
1652 137 : IPosition gridShape(4, nx, ny, npol, nchan);
1653 274 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
1654 274 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1655 137 : IPosition stride(4, 1);
1656 274 : IPosition trc(blc+image->shape()-stride);
1657 137 : LCBox gridBox(blc, trc, gridShape);
1658 274 : SubLattice<Complex> gridSub(*arrayLattice, gridBox);
1659 :
1660 : // Do the copy
1661 137 : image->copyData(gridSub);
1662 137 : }
1663 :
1664 : // IMAGER MIGRATION
1665 : // set sumWeight to 1.0
1666 137 : weights = 1.0;
1667 :
1668 137 : return *image;
1669 : }
1670 :
1671 : // Return weights image
1672 137 : void SDGrid::getWeightImage(ImageInterface<Float>& weightImage, Matrix<Float>& weights)
1673 : {
1674 137 : AlwaysAssert(lattice, AipsError);
1675 137 : AlwaysAssert(image, AipsError);
1676 :
1677 137 : logIO() << LogOrigin("SDGrid", "getWeightImage") << LogIO::NORMAL;
1678 :
1679 137 : weights.resize(sumWeight.shape());
1680 : // IMAGER MIGRATION
1681 : // set sumWeight to 1.0
1682 137 : weights = 1.0;
1683 : // convertArray(weights,sumWeight);
1684 :
1685 137 : weightImage.copyData(*wArrayLattice);
1686 137 : }
1687 :
1688 274 : void SDGrid::ok() {
1689 274 : AlwaysAssert(image, AipsError);
1690 274 : }
1691 :
1692 : // Get the index into the pointing table for this time. Note that the
1693 : // in the pointing table, TIME specifies the beginning of the spanned
1694 : // time range, whereas for the main table, TIME is the centroid.
1695 : // Note that the behavior for multiple matches is not specified! i.e.
1696 : // if there are multiple matches, the index returned depends on the
1697 : // history of previous matches. It is deterministic but not obvious.
1698 : // One could cure this by searching but it would be considerably
1699 : // costlier.
1700 176529 : Int SDGrid::getIndex(const MSPointingColumns& mspc, const Double& time,
1701 : const Double& interval, const Int& antid) {
1702 : //Int start=lastIndex_p;
1703 176529 : Int start=lastIndexPerAnt_p[antid];
1704 176529 : Double tol=interval < 0.0 ? time*C::dbl_epsilon : interval/2.0;
1705 : // Search forwards
1706 176529 : Int nrows=mspc.time().nrow();
1707 957005 : for (Int i=start;i<nrows;i++) {
1708 : // Check for ANTENNA_ID. When antid<0 (default) ANTENNA_ID in POINTING table < 0 is ignored.
1709 : // MS v2 definition indicates ANTENNA_ID in POINING >= 0.
1710 957005 : if (antid >= 0 && mspc.antennaId()(i) != antid) {
1711 293484 : continue;
1712 : }
1713 :
1714 663521 : Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
1715 : // If the interval in the pointing table is negative, use the last
1716 : // entry. Note that this may be invalid (-1) but in that case
1717 : // the calling routine will generate an error
1718 663521 : Double mspc_interval = mspc.interval()(i);
1719 :
1720 663521 : if(mspc_interval<0.0) {
1721 : //return lastIndex_p;
1722 0 : return lastIndexPerAnt_p[antid];
1723 : }
1724 : // Pointing table interval is specified so we have to do a match
1725 : else {
1726 : // Is the midpoint of this pointing table entry within the specified
1727 : // tolerance of the main table entry?
1728 663521 : if(abs(midpoint-time) <= (mspc_interval/2.0+tol)) {
1729 : //lastIndex_p=i;
1730 176529 : lastIndexPerAnt_p[antid]=i;
1731 176529 : return i;
1732 : }
1733 : }
1734 : }
1735 : // Search backwards
1736 0 : for (Int i=start;i>=0;i--) {
1737 0 : if (antid >= 0 && mspc.antennaId()(i) != antid) {
1738 0 : continue;
1739 : }
1740 0 : Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
1741 0 : Double mspc_interval = mspc.interval()(i);
1742 0 : if(mspc_interval<0.0) {
1743 : //return lastIndex_p;
1744 0 : return lastIndexPerAnt_p[antid];
1745 : }
1746 : // Pointing table interval is specified so we have to do a match
1747 : else {
1748 : // Is the midpoint of this pointing table entry within the specified
1749 : // tolerance of the main table entry?
1750 0 : if (abs(midpoint-time) <= (mspc_interval/2.0+tol)) {
1751 : //lastIndex_p=i;
1752 0 : lastIndexPerAnt_p[antid]=i;
1753 0 : return i;
1754 : }
1755 : }
1756 : }
1757 : // No match!
1758 0 : return -1;
1759 : }
1760 :
1761 176511 : Bool SDGrid::getXYPos(const vi::VisBuffer2& vb, Int row) {
1762 :
1763 : // Select the POINTING table (and columns) we'll work with
1764 176511 : const auto haveConvertedColumn = ramPointingTable.nrow() > 0;
1765 :
1766 176511 : const auto & pointingColumns = haveConvertedColumn ?
1767 0 : *ramPointingColumnsPtr
1768 176511 : : vb.subtableColumns().pointing();
1769 :
1770 176511 : const auto nPointings = pointingColumns.nrow();
1771 176511 : const auto havePointings = nPointings > 0;
1772 :
1773 : // We'll need to call these many times, so let's call them once for good
1774 176511 : const auto rowTime = vb.time()(row);
1775 176511 : const auto rowTimeInterval = vb.timeInterval()(row);
1776 176511 : const auto rowAntenna1 = vb.antenna1()(row);
1777 :
1778 : // 1. Try to find the index of a pointing recorded:
1779 : // - for the antenna of the specified row,
1780 : // - at a time close enough to the time at which data was taken
1781 176511 : Int pointingIndex = -1;
1782 :
1783 176511 : if (havePointings) {
1784 : // if (vb.newMS() vb.newMS does not work well using msid
1785 : // Note about above comment:
1786 : // - vb.newMS probably works well
1787 : // - but if the calling code is iterating over the rows of a subchunk
1788 : // vb.newMS returns true for all rows belonging to the first subchunk
1789 : // of the first chunk of a new MS.
1790 : // ???
1791 : // What if vb changed since we were last called ?
1792 : // What if the calling code calls put and get, with different VisBuffers ?
1793 176511 : if (vb.msId() != msId_p) {
1794 155 : lastIndex_p = 0; // No longer used ?
1795 155 : if (lastIndexPerAnt_p.nelements() < (size_t)vb.nAntennas()) {
1796 115 : lastIndexPerAnt_p.resize(vb.nAntennas());
1797 : }
1798 155 : lastIndexPerAnt_p = 0;
1799 155 : msId_p = vb.msId();
1800 155 : lastAntID_p = -1;
1801 : }
1802 : // Try to locate a pointing verifying:
1803 : // | POINTING.TIME - MAIN.TIME | <= 0.5*(POINTING.INTERVAL + tolerance)
1804 : // using first a tiny tolerance, then MAIN.INTERVAL
1805 176511 : constexpr Double useTinyTolerance = -1.0;
1806 176511 : Bool foundPointing {False};
1807 176511 : for(const auto tolerance : {useTinyTolerance, rowTimeInterval}) {
1808 176511 : pointingIndex = getIndex(
1809 : pointingColumns, rowTime, tolerance , rowAntenna1
1810 : );
1811 176511 : foundPointing = pointingIndex >= 0;
1812 176511 : if (foundPointing) break;
1813 : }
1814 :
1815 : // Making the implicit type conversion explicit.
1816 : // Conversion is safe because it occurs only when pointingIndex >= 0.
1817 176511 : const auto foundValidPointing = (
1818 176511 : foundPointing and (static_cast<rownr_t>(pointingIndex) < nPointings)
1819 : );
1820 :
1821 176511 : if (not foundValidPointing) {
1822 0 : LogIO logger(LogOrigin("SDGrid","getXYPos"));
1823 0 : logger << LogIO::DEBUGGING;
1824 0 : logger.output()
1825 0 : << "Failed to find pointing information for time "
1826 0 : << MVTime(rowTime/86400.0)
1827 0 : << " : omitting this point";
1828 0 : logger << LogIO::POST;
1829 0 : return false;
1830 0 : }
1831 : }
1832 :
1833 : // 2. At this stage we have:
1834 : // * either no pointings
1835 : // * or pointings and a valid pointingIndex
1836 : // Decide now if we need to interpolate antenna's pointing direction
1837 : // at data-taking time:
1838 : // we'll do so when data is sampled faster than pointings are recorded
1839 176511 : Bool needInterpolation = False;
1840 176511 : if (havePointings) {
1841 176511 : const auto pointingInterval = pointingColumns.interval()(pointingIndex);
1842 176511 : if (rowTimeInterval < pointingInterval) needInterpolation = True;
1843 : }
1844 176511 : const auto mustInterpolate = havePointings && needInterpolation;
1845 :
1846 : // 3. Create interpolator if needed
1847 176511 : if (mustInterpolate) {
1848 19602 : if (not isSplineInterpolationReady) {
1849 : const auto nAntennas = static_cast<size_t>(
1850 16 : vb.ms().antenna().nrow()
1851 : );
1852 16 : interpolator = new SDPosInterpolator(
1853 : pointingColumns,
1854 16 : pointingDirCol_p,
1855 : nAntennas
1856 16 : );
1857 16 : isSplineInterpolationReady = true;
1858 : } else {
1859 : // We have an interpolator. Re-use it if possible.
1860 : const auto canReuseInterpolator =
1861 19586 : interpolator->inTimeRange(rowTime, rowAntenna1);
1862 19586 : if (not canReuseInterpolator) {
1863 : // setup spline interpolator for the current dataset
1864 : // (CAS-11261, 2018/5/22 WK)
1865 : // delete and re-create it
1866 188 : delete interpolator;
1867 188 : interpolator = 0;
1868 : const auto nAntennas = static_cast<size_t>(
1869 188 : vb.ms().antenna().nrow()
1870 : );
1871 188 : interpolator = new SDPosInterpolator(
1872 : pointingColumns,
1873 188 : pointingDirCol_p,
1874 : nAntennas
1875 188 : );
1876 : }
1877 : }
1878 : }
1879 :
1880 : // 4. Create the direction conversion machine if needed
1881 353022 : if ( pointingDirCol_p == "SOURCE_OFFSET" or
1882 176511 : pointingDirCol_p == "POINTING_OFFSET" ) {
1883 : // it makes no sense to track in offset coordinates...
1884 : // hopefully the user set the image coords right
1885 0 : fixMovingSource_p = false;
1886 : }
1887 :
1888 176511 : const auto needDirectionConverter = (
1889 176511 : not havePointings or not haveConvertedColumn or fixMovingSource_p
1890 : );
1891 :
1892 176511 : if (not pointingToImage and needDirectionConverter) {
1893 : // Setup our Measures container
1894 : const auto & rowAntenna1Position =
1895 132 : vb.subtableColumns().antenna().positionMeas()(rowAntenna1);
1896 : // Set dummy time stamp 1 day before rowTime
1897 264 : const MEpoch dummyEpoch(Quantity(rowTime - 86400.0, "s"));
1898 : // mFrame_p = MeasFrame(dummyEpoch, rowAntenna1Position);
1899 132 : mFrame_p.epoch() ? mFrame_p.resetEpoch(dummyEpoch)
1900 0 : : mFrame_p.set(dummyEpoch);
1901 132 : mFrame_p.position() ? mFrame_p.resetPosition(rowAntenna1Position)
1902 0 : : mFrame_p.set(rowAntenna1Position);
1903 : // Remember antenna id for next call,
1904 : // which may be done using a different VisBuffer ...
1905 132 : lastAntID_p = rowAntenna1;
1906 : // Compute the "model" required to setup the direction conversion machine
1907 132 : if (havePointings) {
1908 264 : worldPosMeas = mustInterpolate ?
1909 : directionMeas(pointingColumns, pointingIndex, rowTime)
1910 132 : : directionMeas(pointingColumns, pointingIndex);
1911 : } else {
1912 : // Without pointings, this sets the direction to the phase center
1913 0 : worldPosMeas = vb.direction1()(row);
1914 : }
1915 : // Make a direction conversion machine, converting
1916 : // from: the reference frame of the "model"
1917 : // to: image's reference frame
1918 132 : MDirection::Ref outRef(directionCoord.directionType(), mFrame_p);
1919 132 : pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
1920 132 : if (not pointingToImage) {
1921 0 : logIO_p << "Cannot make direction conversion machine" << LogIO::EXCEPTION;
1922 : }
1923 : // Perform 1 dummy direction conversion to clear values
1924 : // cached in static variables of casacore functions like MeasTable::dUT1
1925 132 : MDirection _dir_tmp = (*pointingToImage)();
1926 132 : }
1927 :
1928 353022 : const MEpoch rowEpoch(Quantity(rowTime, "s"));
1929 : { // 5. Update the frame holding the measurements for this row
1930 : // ---- Always reset the epoch
1931 176511 : mFrame_p.resetEpoch(rowEpoch);
1932 : // ---- Reset antenna position only if antenna changed
1933 : // since we were last called
1934 176511 : const auto antennaChanged = (lastAntID_p != rowAntenna1);
1935 176511 : if (antennaChanged) {
1936 : { // Debug messages
1937 14588 : if (lastAntID_p == -1) {
1938 : // antenna ID is unset
1939 23 : logIO_p << LogIO::DEBUGGING
1940 : << "updating antenna position for conversion: new MS ID " << msId_p
1941 23 : << ", antenna ID " << rowAntenna1 << LogIO::POST;
1942 : } else {
1943 14565 : logIO_p << LogIO::DEBUGGING
1944 : << "updating antenna position for conversion: MS ID " << msId_p
1945 : << ", last antenna ID " << lastAntID_p
1946 14565 : << ", new antenna ID " << rowAntenna1 << LogIO::POST;
1947 : }
1948 : }
1949 : MPosition rowAntenna1Position (
1950 14588 : vb.subtableColumns().antenna().positionMeas()(rowAntenna1)
1951 14588 : );
1952 14588 : mFrame_p.resetPosition(rowAntenna1Position);
1953 : // Remember antenna id for next call,
1954 : // which may be done using a different VisBuffer ...
1955 14588 : lastAntID_p = rowAntenna1;
1956 14588 : }
1957 : }
1958 :
1959 : // 6. Compute user-specified column direction at data-taking time,
1960 : // converted to image's direction reference frame
1961 176511 : if (havePointings) {
1962 : const auto columnDirection = mustInterpolate ?
1963 : directionMeas(pointingColumns, pointingIndex, rowTime)
1964 176511 : : directionMeas(pointingColumns, pointingIndex);
1965 : worldPosMeas = haveConvertedColumn ?
1966 : columnDirection
1967 176511 : : (*pointingToImage)(columnDirection);
1968 : { // Old debug stuff
1969 : //Vector<Double> newdirv = newdir.getAngle("rad").getValue();
1970 : //cerr<<"dir0="<<newdirv(0)<<endl;
1971 : //fprintf(pfile,"%.8f %.8f \n", newdirv(0), newdirv(1));
1972 : //printf("%lf %lf \n", newdirv(0), newdirv(1));
1973 : }
1974 176511 : } else {
1975 : // Without pointings, this converts the direction of the phase center ?
1976 0 : worldPosMeas = (*pointingToImage)(vb.direction1()(row));
1977 : }
1978 :
1979 : // 7. Convert world direction coordinates to image pixel coordinates
1980 176511 : Bool havePixel = directionCoord.toPixel(xyPos, worldPosMeas);
1981 176511 : if (not havePixel) { // Log warning
1982 0 : logIO_p << LogIO::WARN
1983 : << "Failed to find a pixel for pointing direction of "
1984 0 : << MVTime(worldPosMeas.getValue().getLong("rad")).string(MVTime::TIME)
1985 : << ", "
1986 0 : << MVAngle(worldPosMeas.getValue().getLat("rad")).string(MVAngle::ANGLE)
1987 0 : << LogIO::POST;
1988 0 : return false;
1989 : }
1990 :
1991 : // 8. Handle moving sources
1992 176511 : if (fixMovingSource_p) {
1993 21292 : if (xyPosMovingOrig_p.nelements() < 2) {
1994 5 : directionCoord.toPixel(xyPosMovingOrig_p, firstMovingDir_p);
1995 : }
1996 : //via HADEC or AZEL for parallax of near sources
1997 21292 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
1998 21292 : MDirection tmphadec;
1999 21292 : if (upcase(movingDir_p.getRefString()).contains("APP")) {
2000 11656 : tmphadec = MDirection::Convert(
2001 11656 : vbutil_p->getEphemDir(vb, phaseCenterTime_p),
2002 : outref1
2003 5828 : )();
2004 15464 : } else if (upcase(movingDir_p.getRefString()).contains("COMET")) {
2005 2914 : MeasComet mcomet(Path(ephemTableName_p).absoluteName());
2006 2914 : mFrame_p.set(mcomet);
2007 2914 : tmphadec = MDirection::Convert(movingDir_p, outref1)();
2008 2914 : } else {
2009 12550 : tmphadec = MDirection::Convert(movingDir_p, outref1)();
2010 : }
2011 21292 : MDirection actSourceDir = (*pointingToImage)(tmphadec);
2012 21292 : Vector<Double> actPix;
2013 21292 : directionCoord.toPixel(actPix, actSourceDir);
2014 :
2015 : // cout << row
2016 : // << " scan " << vb.scan()(row)
2017 : // << " xyPos " << xyPos
2018 : // << " xyposmovorig " << xyPosMovingOrig_p
2019 : // << " actPix " << actPix << endl;
2020 :
2021 21292 : xyPos = xyPos + xyPosMovingOrig_p - actPix;
2022 21292 : }
2023 :
2024 176511 : return havePixel;
2025 176511 : }
2026 :
2027 157061 : MDirection SDGrid::directionMeas(const MSPointingColumns& mspc, const Int& index){
2028 157061 : if (pointingDirCol_p == "TARGET") {
2029 0 : return mspc.targetMeas(index);
2030 157061 : } else if (pointingDirCol_p == "POINTING_OFFSET") {
2031 0 : if (!mspc.pointingOffsetMeasCol().isNull()) {
2032 0 : return mspc.pointingOffsetMeas(index);
2033 : }
2034 0 : cerr << "No PONTING_OFFSET column in POINTING table" << endl;
2035 157061 : } else if (pointingDirCol_p == "SOURCE_OFFSET") {
2036 0 : if (!mspc.sourceOffsetMeasCol().isNull()) {
2037 0 : return mspc.sourceOffsetMeas(index);
2038 : }
2039 0 : cerr << "No SOURCE_OFFSET column in POINTING table" << endl;
2040 157061 : } else if (pointingDirCol_p == "ENCODER") {
2041 0 : if (!mspc.encoderMeas().isNull()) {
2042 0 : return mspc.encoderMeas()(index);
2043 : }
2044 0 : cerr << "No ENCODER column in POINTING table" << endl;
2045 : }
2046 :
2047 : //default return this
2048 157061 : return mspc.directionMeas(index);
2049 : }
2050 :
2051 : // for the cases, interpolation of the pointing direction requires
2052 : // when data sampling rate higher than the pointing data recording
2053 : // (e.g. fast OTF)
2054 19618 : MDirection SDGrid::directionMeas(const MSPointingColumns& mspc, const Int& index,
2055 : const Double& time){
2056 : //spline interpolation
2057 19618 : if (isSplineInterpolationReady) {
2058 19618 : Int antid = mspc.antennaId()(index);
2059 19618 : if (interpolator->doSplineInterpolation(antid)) {
2060 19618 : return interpolator->interpolateDirectionMeasSpline(mspc, time, index, antid);
2061 : }
2062 : }
2063 :
2064 : //linear interpolation (as done before CAS-7911)
2065 : Int index1, index2;
2066 0 : if (time < mspc.time()(index)) {
2067 0 : if (index > 0) {
2068 0 : index1 = index-1;
2069 0 : index2 = index;
2070 0 : } else if (index == 0) {
2071 0 : index1 = index;
2072 0 : index2 = index+1;
2073 : }
2074 : } else {
2075 0 : if (index < Int(mspc.nrow()-1)) {
2076 0 : index1 = index;
2077 0 : index2 = index+1;
2078 0 : } else if (index == Int(mspc.nrow()-1) || (mspc.time()(index)-mspc.time()(index+1))>2*mspc.interval()(index)) {
2079 0 : index1 = index-1;
2080 0 : index2 = index;
2081 : }
2082 : }
2083 0 : return interpolateDirectionMeas(mspc, time, index, index1, index2);
2084 : }
2085 :
2086 0 : MDirection SDGrid::interpolateDirectionMeas(const MSPointingColumns& mspc,
2087 : const Double& time,
2088 : const Int& index,
2089 : const Int& indx1,
2090 : const Int& indx2){
2091 0 : Vector<Double> dir1,dir2;
2092 0 : Vector<Double> newdir(2),scanRate(2);
2093 : Double dLon, dLat;
2094 : Double ftime,ftime2,ftime1,dtime;
2095 0 : MDirection newDirMeas;
2096 0 : MDirection::Ref rf;
2097 : Bool isfirstRefPt;
2098 :
2099 0 : if (indx1 == index) {
2100 0 : isfirstRefPt = true;
2101 : } else {
2102 0 : isfirstRefPt = false;
2103 : }
2104 :
2105 0 : if (pointingDirCol_p == "TARGET") {
2106 0 : dir1 = mspc.targetMeas(indx1).getAngle("rad").getValue();
2107 0 : dir2 = mspc.targetMeas(indx2).getAngle("rad").getValue();
2108 0 : } else if (pointingDirCol_p == "POINTING_OFFSET") {
2109 0 : if (!mspc.pointingOffsetMeasCol().isNull()) {
2110 0 : dir1 = mspc.pointingOffsetMeas(indx1).getAngle("rad").getValue();
2111 0 : dir2 = mspc.pointingOffsetMeas(indx2).getAngle("rad").getValue();
2112 : } else {
2113 0 : cerr << "No PONTING_OFFSET column in POINTING table" << endl;
2114 : }
2115 0 : } else if (pointingDirCol_p == "SOURCE_OFFSET") {
2116 0 : if (!mspc.sourceOffsetMeasCol().isNull()) {
2117 0 : dir1 = mspc.sourceOffsetMeas(indx1).getAngle("rad").getValue();
2118 0 : dir2 = mspc.sourceOffsetMeas(indx2).getAngle("rad").getValue();
2119 : } else {
2120 0 : cerr << "No SOURCE_OFFSET column in POINTING table" << endl;
2121 : }
2122 0 : } else if (pointingDirCol_p == "ENCODER") {
2123 0 : if (!mspc.encoderMeas().isNull()) {
2124 0 : dir1 = mspc.encoderMeas()(indx1).getAngle("rad").getValue();
2125 0 : dir2 = mspc.encoderMeas()(indx2).getAngle("rad").getValue();
2126 : } else {
2127 0 : cerr << "No ENCODER column in POINTING table" << endl;
2128 : }
2129 : } else {
2130 0 : dir1 = mspc.directionMeas(indx1).getAngle("rad").getValue();
2131 0 : dir2 = mspc.directionMeas(indx2).getAngle("rad").getValue();
2132 : }
2133 :
2134 0 : dLon = dir2(0) - dir1(0);
2135 0 : dLat = dir2(1) - dir1(1);
2136 0 : ftime = floor(mspc.time()(indx1));
2137 0 : ftime2 = mspc.time()(indx2) - ftime;
2138 0 : ftime1 = mspc.time()(indx1) - ftime;
2139 0 : dtime = ftime2 - ftime1;
2140 0 : scanRate(0) = dLon/dtime;
2141 0 : scanRate(1) = dLat/dtime;
2142 : //scanRate(0) = dir2(0)/dtime-dir1(0)/dtime;
2143 : //scanRate(1) = dir2(1)/dtime-dir1(1)/dtime;
2144 : //Double delT = mspc.time()(index)-time;
2145 : //cerr<<"index="<<index<<" dLat="<<dLat<<" dtime="<<dtime<<" delT="<< delT<<endl;
2146 : //cerr<<"deldirlat="<<scanRate(1)*fabs(delT)<<endl;
2147 0 : if (isfirstRefPt) {
2148 0 : newdir(0) = dir1(0)+scanRate(0)*fabs(mspc.time()(index)-time);
2149 0 : newdir(1) = dir1(1)+scanRate(1)*fabs(mspc.time()(index)-time);
2150 0 : rf = mspc.directionMeas(indx1).getRef();
2151 : } else {
2152 0 : newdir(0) = dir2(0)-scanRate(0)*fabs(mspc.time()(index)-time);
2153 0 : newdir(1) = dir2(1)-scanRate(1)*fabs(mspc.time()(index)-time);
2154 0 : rf = mspc.directionMeas(indx2).getRef();
2155 : }
2156 : //default return this
2157 0 : Quantity rDirLon(newdir(0), "rad");
2158 0 : Quantity rDirLat(newdir(1), "rad");
2159 0 : newDirMeas = MDirection(rDirLon, rDirLat, rf);
2160 : //cerr<<"newDirMeas rf="<<newDirMeas.getRefString()<<endl;
2161 : //return mspc.directionMeas(index);
2162 0 : return newDirMeas;
2163 0 : }
2164 :
2165 169226 : void SDGrid::pickWeights(const vi::VisBuffer2& vb, Matrix<Float>& weight){
2166 : //break reference
2167 169226 : weight.resize();
2168 :
2169 169226 : if (useImagingWeight_p) {
2170 0 : weight.reference(vb.imagingWeight());
2171 : } else {
2172 169226 : const Cube<Float> weightSpec(vb.weightSpectrum());
2173 169226 : weight.resize(vb.nChannels(), vb.nRows());
2174 :
2175 : // CAS-9957 correct weight propagation from linear/circular correlations to Stokes I
2176 94633912 : const auto toStokesWeight = [](float weight0, float weight1) {
2177 94633912 : const auto denominator = weight0 + weight1;
2178 94633912 : const auto numerator = weight0 * weight1;
2179 94633912 : constexpr float fmin = std::numeric_limits<float>::min();
2180 94633912 : return abs(denominator) < fmin ? 0.0f : 4.0f * numerator / denominator;
2181 : };
2182 :
2183 169226 : if (weightSpec.nelements() == 0) {
2184 0 : const auto &weightMat = vb.weight();
2185 0 : const ssize_t npol = weightMat.shape()(0);
2186 0 : if (npol == 1) {
2187 0 : const auto weight0 = weightMat.row(0);
2188 0 : for (rownr_t k = 0; k < vb.nRows(); ++k) {
2189 0 : weight.column(k).set(weight0(k));
2190 : }
2191 0 : } else if (npol == 2) {
2192 0 : const auto weight0 = weightMat.row(0);
2193 0 : const auto weight1 = weightMat.row(1);
2194 0 : for (rownr_t k = 0; k < vb.nRows(); ++k) {
2195 : //cerr << "nrow " << vb.nRow() << " " << weight.shape() << " " << weight.column(k).shape() << endl;
2196 0 : weight.column(k).set(toStokesWeight(weight0(k), weight1(k)));
2197 : }
2198 0 : } else {
2199 : // It seems current code doesn't support 4 pol case. So, give up
2200 : // processing such data to avoid producing unintended result
2201 0 : throw AipsError("Imaging full-Stokes data (npol=4) is not supported.");
2202 : }
2203 : } else {
2204 169226 : const ssize_t npol = weightSpec.shape()(0);
2205 169226 : if (npol == 1) {
2206 50 : weight = weightSpec.yzPlane(0);
2207 169176 : } else if (npol == 2) {
2208 169176 : const auto weight0 = weightSpec.yzPlane(0);
2209 169176 : const auto weight1 = weightSpec.yzPlane(1);
2210 345637 : for (rownr_t k = 0; k < vb.nRows(); ++k) {
2211 94810373 : for (int chan = 0; chan < vb.nChannels(); ++chan) {
2212 94633912 : weight(chan, k) = toStokesWeight(weight0(chan, k), weight1(chan, k));
2213 : }
2214 : }
2215 169176 : } else {
2216 : // It seems current code doesn't support 4 pol case. So, give up
2217 : // processing such data to avoid producing unintended result
2218 0 : throw AipsError("Imaging full-Stokes data (npol=4) is not supported.");
2219 : }
2220 : }
2221 169226 : }
2222 169226 : }
2223 :
2224 137 : void SDGrid::clipMinMax() {
2225 137 : if (clipminmax_) {
2226 : Bool gmin_b, gmax_b, wmin_b, wmax_b, np_b;
2227 8 : const auto *gmin_p = gmin_.getStorage(gmin_b);
2228 8 : const auto *gmax_p = gmax_.getStorage(gmax_b);
2229 8 : const auto *wmin_p = wmin_.getStorage(wmin_b);
2230 8 : const auto *wmax_p = wmax_.getStorage(wmax_b);
2231 8 : const auto *np_p = npoints_.getStorage(np_b);
2232 :
2233 : Bool data_b, weight_b, sumw_b;
2234 8 : auto data_p = griddedData.getStorage(data_b);
2235 8 : auto weight_p = wGriddedData.getStorage(weight_b);
2236 8 : auto sumw_p = sumWeight.getStorage(sumw_b);
2237 :
2238 8 : auto arrayShape = griddedData.shape();
2239 8 : size_t num_xy = arrayShape.getFirst(2).product();
2240 8 : size_t num_polchan = arrayShape.getLast(2).product();
2241 80 : for (size_t i = 0; i < num_xy; ++i) {
2242 153 : for (size_t j = 0; j < num_polchan; ++j) {
2243 81 : auto k = i * num_polchan + j;
2244 81 : if (np_p[k] == 1) {
2245 2 : auto wt = wmin_p[k];
2246 2 : data_p[k] = wt * gmin_p[k];
2247 2 : weight_p[k] = wt;
2248 2 : sumw_p[j] += wt;
2249 79 : } else if (np_p[k] == 2) {
2250 1 : auto wt = wmin_p[k];
2251 1 : data_p[k] = wt * gmin_p[k];
2252 1 : weight_p[k] = wt;
2253 1 : sumw_p[j] += wt;
2254 1 : wt = wmax_p[k];
2255 1 : data_p[k] += wt * gmax_p[k];
2256 1 : weight_p[k] += wt;
2257 1 : sumw_p[j] += wt;
2258 : }
2259 : }
2260 : }
2261 :
2262 8 : wGriddedData.putStorage(weight_p, weight_b);
2263 8 : griddedData.putStorage(data_p, data_b);
2264 8 : sumWeight.putStorage(sumw_p, sumw_b);
2265 :
2266 8 : npoints_.freeStorage(np_p, np_b);
2267 8 : wmax_.freeStorage(wmax_p, wmax_b);
2268 8 : wmin_.freeStorage(wmin_p, wmin_b);
2269 8 : gmax_.freeStorage(gmax_p, gmax_b);
2270 8 : gmin_.freeStorage(gmin_p, gmin_b);
2271 8 : }
2272 137 : }
2273 :
2274 274 : const String & SDGrid::toString(const ConvertFirst convertFirst) {
2275 : static const std::array<String,3> name {
2276 : "never",
2277 : "always",
2278 : "auto"
2279 274 : };
2280 :
2281 274 : switch (convertFirst) {
2282 274 : case ConvertFirst::NEVER:
2283 : case ConvertFirst::ALWAYS:
2284 : case ConvertFirst::AUTO:
2285 274 : return name[static_cast<size_t>(convertFirst)];
2286 0 : default:
2287 0 : String errMsg {"Illegal ConvertFirst enum: "};
2288 0 : errMsg += String::toString(static_cast<Int>(convertFirst));
2289 0 : throw AipsError(
2290 : errMsg,
2291 : __FILE__,
2292 : __LINE__,
2293 : AipsError::Category::INVALID_ARGUMENT
2294 0 : );
2295 : // Avoid potential compiler warning
2296 : return name[static_cast<size_t>(ConvertFirst::NEVER)];
2297 : }
2298 : }
2299 :
2300 274 : SDGrid::ConvertFirst SDGrid::fromString(const String & name) {
2301 : static const std::array<ConvertFirst,3> schemes {
2302 : ConvertFirst::NEVER,
2303 : ConvertFirst::ALWAYS,
2304 : ConvertFirst::AUTO
2305 : };
2306 :
2307 274 : for (const auto scheme : schemes) {
2308 274 : if (name == toString(scheme)) return scheme;
2309 : }
2310 :
2311 0 : String errMsg {"Illegal ConvertFirst name: "};
2312 0 : errMsg += name;
2313 0 : throw AipsError(
2314 : errMsg,
2315 : __FILE__,
2316 : __LINE__,
2317 : AipsError::Category::INVALID_ARGUMENT
2318 0 : );
2319 : // Avoid potential compiler warning
2320 : return ConvertFirst::NEVER;
2321 0 : }
2322 :
2323 274 : void SDGrid::setConvertFirst(const String &name) {
2324 274 : convertFirst = fromString(name);
2325 274 : }
2326 :
2327 : } // End of namespace: refim
2328 : } // End of namespace: casa
|
