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: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/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/VisBuffer.h>
75 : #include <msvis/MSVis/VisibilityIterator.h>
76 : #include <casacore/scimath/Mathematics/RigidVector.h>
77 : #include <synthesis/MeasurementComponents/SDGrid.h>
78 : #include <synthesis/TransformMachines/SkyJones.h>
79 : #include <synthesis/TransformMachines/StokesImageUtil.h>
80 :
81 : #if defined(SDGRID_PERFS)
82 : #include <iomanip>
83 : #include <exception>
84 : #include <sstream>
85 : #endif
86 :
87 : using namespace casacore;
88 : namespace casa {
89 :
90 : #if defined(SDGRID_PERFS)
91 : namespace sdgrid_perfs {
92 1612 : ChronoStat::ChronoStat(const string & name)
93 : : name_ {name},
94 : started_ {false},
95 : n_laps_ {0},
96 : n_overflows_ {0},
97 : n_underflows_ {0},
98 1612 : laps_sum_ {Duration::zero()},
99 1612 : laps_min_ {Duration::max()},
100 3224 : laps_max_ {Duration::min()}
101 1612 : {}
102 :
103 1612 : const std::string& ChronoStat::name() const {
104 1612 : return name_;
105 : }
106 :
107 1612 : void ChronoStat::setName(const std::string& name) {
108 1612 : name_ = name;
109 1612 : }
110 :
111 571954 : void ChronoStat::start() {
112 571954 : if (not started_) {
113 571954 : started_ = true;
114 571954 : ++n_laps_;
115 571954 : lap_start_time_ = Clock::now();
116 : } else {
117 0 : ++n_overflows_;
118 : }
119 571954 : }
120 571954 : void ChronoStat::stop() {
121 571954 : if (started_) {
122 571954 : auto lap_duration = Clock::now() - lap_start_time_;
123 571954 : laps_sum_ += lap_duration;
124 571954 : started_ = false;
125 571954 : if (lap_duration < laps_min_) laps_min_ = lap_duration;
126 571954 : if (lap_duration > laps_max_) laps_max_ = lap_duration;
127 : } else {
128 0 : ++n_underflows_;
129 : }
130 571954 : }
131 :
132 1612 : ChronoStat::Duration ChronoStat::lapsSum() const {
133 1612 : return laps_sum_;
134 : }
135 :
136 1612 : ChronoStat::Duration ChronoStat::lapsMin() const {
137 1612 : return n_laps_ > 0 ? laps_min_ : Duration::zero();
138 : }
139 :
140 1612 : ChronoStat::Duration ChronoStat::lapsMax() const {
141 1612 : return n_laps_ > 0 ? laps_max_ : Duration::zero();
142 : }
143 :
144 1612 : ChronoStat::Duration ChronoStat::lapsMean() const {
145 1612 : return n_laps_ > 0 ? laps_sum_ / n_laps_ : Duration::zero();
146 : }
147 :
148 1612 : unsigned int ChronoStat::lapsCount() const {
149 1612 : return n_laps_;
150 : }
151 :
152 0 : bool ChronoStat::isEmpty() const {
153 0 : return n_laps_ == 0;
154 : }
155 :
156 1612 : unsigned int ChronoStat::nOverflows() const {
157 1612 : return n_overflows_;
158 : }
159 :
160 1612 : unsigned int ChronoStat::nUnderflows() const {
161 1612 : return n_underflows_;
162 : }
163 :
164 12896 : std::string ChronoStat::quote(const std::string& s) const {
165 12896 : return "\"" + s + "\"";
166 : }
167 :
168 1612 : std::string ChronoStat::json() const {
169 3224 : std::ostringstream os;
170 1612 : os << quote(name()) << ": {"
171 3224 : << quote("sum") << ": " << lapsSum().count()
172 3224 : << " ," << quote("count") << ": " << lapsCount()
173 3224 : << " ," << quote("min") << ": " << lapsMin().count()
174 3224 : << " ," << quote("mean") << ": " << lapsMean().count()
175 3224 : << " ," << quote("max") << ": " << lapsMax().count()
176 3224 : << " ," << quote("overflows") << ": " << nOverflows()
177 3224 : << " ," << quote("underflows") << ": " << nUnderflows()
178 1612 : << "}";
179 3224 : return os.str();
180 : }
181 :
182 0 : std::ostream& operator<<(std::ostream &os, const ChronoStat &c) {
183 0 : constexpr auto eol = '\n';
184 0 : os << "name: " << c.name() << eol
185 0 : << "sum: " << std::setw(20) << std::right << c.lapsSum().count() << eol
186 0 : << "count: " << std::setw(20) << std::right << c.lapsCount() << eol
187 0 : << "min: " << std::setw(20) << std::right << c.lapsMin().count() << eol
188 0 : << "mean: " << std::setw(20) << std::right << c.lapsMean().count() << eol
189 0 : << "max: " << std::setw(20) << std::right << c.lapsMax().count() << eol
190 0 : << "overflows: " << std::setw(20) << std::right << c.nOverflows() << eol
191 0 : << "underflows: " << std::setw(20) << std::right << c.nUnderflows();
192 0 : return os;
193 : }
194 :
195 561954 : StartStop::StartStop(ChronoStat &c)
196 561954 : : c_ {c} {
197 561954 : c_.start();
198 561954 : }
199 :
200 1123908 : StartStop::~StartStop() {
201 561954 : c_.stop();
202 561954 : }
203 :
204 : } // namespace sdgrid_perfs
205 :
206 : using namespace sdgrid_perfs;
207 :
208 : #endif
209 :
210 10 : SDGrid::SDGrid(SkyJones& sj, Int icachesize, Int itilesize,
211 10 : String iconvType, Int userSupport, Bool useImagingWeight)
212 : : FTMachine(), sj_p(&sj), imageCache(0), wImageCache(0),
213 : cachesize(icachesize), tilesize(itilesize),
214 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
215 10 : pointingToImage(0), rowPixel {xyPos, false}, userSetSupport_p(userSupport),
216 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
217 : minWeight_p(0.), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
218 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(false),
219 : cache {Cache(*(const_cast<SDGrid *>(this)))},
220 10 : cacheIsEnabled {false}
221 : {
222 10 : lastIndex_p=0;
223 10 : initPerfs();
224 10 : }
225 :
226 18 : SDGrid::SDGrid(MPosition& mLocation, SkyJones& sj, Int icachesize, Int itilesize,
227 18 : String iconvType, Int userSupport, Float minweight, Bool clipminmax, Bool useImagingWeight)
228 : : FTMachine(), sj_p(&sj), imageCache(0), wImageCache(0),
229 : cachesize(icachesize), tilesize(itilesize),
230 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
231 18 : pointingToImage(0), rowPixel {xyPos, false}, userSetSupport_p(userSupport),
232 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
233 : minWeight_p(minweight), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
234 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(clipminmax),
235 : cache {Cache(*(const_cast<SDGrid *>(this)))},
236 18 : cacheIsEnabled {false}
237 : {
238 18 : mLocation_p=mLocation;
239 18 : lastIndex_p=0;
240 18 : initPerfs();
241 18 : }
242 :
243 0 : SDGrid::SDGrid(Int icachesize, Int itilesize,
244 0 : String iconvType, Int userSupport, Bool useImagingWeight)
245 : : FTMachine(), sj_p(0), imageCache(0), wImageCache(0),
246 : cachesize(icachesize), tilesize(itilesize),
247 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
248 0 : pointingToImage(0), rowPixel {xyPos, false}, userSetSupport_p(userSupport),
249 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
250 : minWeight_p(0.), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
251 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(false),
252 : cache {Cache(*(const_cast<SDGrid *>(this)))},
253 0 : cacheIsEnabled {false}
254 : {
255 0 : lastIndex_p=0;
256 0 : initPerfs();
257 0 : }
258 :
259 93 : SDGrid::SDGrid(MPosition &mLocation, Int icachesize, Int itilesize,
260 93 : String iconvType, Int userSupport, Float minweight, Bool clipminmax, Bool useImagingWeight)
261 : : FTMachine(), sj_p(0), imageCache(0), wImageCache(0),
262 : cachesize(icachesize), tilesize(itilesize),
263 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
264 93 : pointingToImage(0), rowPixel {xyPos, false}, userSetSupport_p(userSupport),
265 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
266 : minWeight_p(minweight), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1),
267 : msId_p(-1),
268 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(clipminmax),
269 : cache {Cache(*(const_cast<SDGrid *>(this)))},
270 93 : cacheIsEnabled {false}
271 : {
272 93 : mLocation_p=mLocation;
273 93 : lastIndex_p=0;
274 93 : initPerfs();
275 93 : }
276 :
277 3 : SDGrid::SDGrid(MPosition &mLocation, Int icachesize, Int itilesize,
278 : String iconvType, Float truncate, Float gwidth, Float jwidth,
279 3 : Float minweight, Bool clipminmax, Bool useImagingWeight)
280 : : FTMachine(), sj_p(0), imageCache(0), wImageCache(0),
281 : cachesize(icachesize), tilesize(itilesize),
282 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
283 3 : pointingToImage(0), rowPixel {xyPos, false}, userSetSupport_p(-1),
284 : truncate_p(truncate), gwidth_p(gwidth), jwidth_p(jwidth),
285 : minWeight_p(minweight), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
286 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(clipminmax),
287 : cache {Cache(*(const_cast<SDGrid *>(this)))},
288 3 : cacheIsEnabled {false}
289 : {
290 3 : mLocation_p=mLocation;
291 3 : lastIndex_p=0;
292 3 : initPerfs();
293 3 : }
294 :
295 : //----------------------------------------------------------------------
296 :
297 124 : void SDGrid::initPerfs() {
298 : #if defined(SDGRID_PERFS)
299 124 : cNextChunk.setName("iterateNextChunk");
300 124 : cMatchAllSpwChans.setName("matchAllSpwChans");
301 124 : cMatchChannel.setName("matchChannel");
302 124 : cPickWeights.setName("pickWeights");
303 124 : cInterpolateFrequencyToGrid.setName("interpolateFrequencyToGrid");
304 124 : cSearchValidPointing.setName("searchValidPointing");
305 124 : cComputeSplines.setName("computeSplines");
306 124 : cResetFrame.setName("resetFrame");
307 124 : cInterpolateDirection.setName("interpolateDirection");
308 124 : cConvertDirection.setName("convertDirection");
309 124 : cComputeDirectionPixel.setName("computeDirectionPixel");
310 124 : cHandleMovingSource.setName("handleMovingSource");
311 124 : cGridData.setName("gridData");
312 : #endif
313 124 : }
314 :
315 :
316 :
317 : //----------------------------------------------------------------------
318 0 : SDGrid& SDGrid::operator=(const SDGrid& other)
319 : {
320 0 : if(this!=&other) {
321 : //Do the base parameters
322 0 : FTMachine::operator=(other);
323 0 : sj_p=other.sj_p;
324 0 : imageCache=other.imageCache;
325 0 : wImage=other.wImage;
326 0 : wImageCache=other.wImageCache;
327 0 : cachesize=other.cachesize;
328 0 : tilesize=other.tilesize;
329 0 : isTiled=other.isTiled;
330 0 : lattice=other.lattice;
331 0 : arrayLattice=other.arrayLattice;
332 0 : wLattice=other.wLattice;
333 0 : wArrayLattice=other.wArrayLattice;
334 0 : convType=other.convType;
335 0 : if(other.wImage !=0)
336 0 : wImage=(TempImage<Float> *)other.wImage->cloneII();
337 : else
338 0 : wImage=0;
339 0 : pointingDirCol_p=other.pointingDirCol_p;
340 0 : pointingToImage=0;
341 0 : xyPos.resize();
342 0 : xyPos=other.xyPos;
343 0 : rowPixel = MaskedPixelRef(xyPos, false);
344 0 : xyPosMovingOrig_p=other.xyPosMovingOrig_p;
345 0 : convFunc.resize();
346 0 : convFunc=other.convFunc;
347 0 : convSampling=other.convSampling;
348 0 : convSize=other.convSize;
349 0 : convSupport=other.convSupport;
350 0 : userSetSupport_p=other.userSetSupport_p;
351 0 : lastIndex_p=0;
352 0 : lastIndexPerAnt_p=0;
353 0 : lastAntID_p=-1;
354 0 : msId_p=-1;
355 0 : useImagingWeight_p=other.useImagingWeight_p;
356 0 : clipminmax_=other.clipminmax_;
357 0 : cacheIsEnabled=false;
358 0 : cache=Cache(*(const_cast<SDGrid *>(this)));
359 : };
360 0 : return *this;
361 : };
362 :
363 474 : String SDGrid::name() const{
364 474 : return String("SDGrid");
365 : }
366 :
367 : void
368 114 : SDGrid::setEnableCache(Bool doEnable) {
369 114 : cacheIsEnabled = doEnable;
370 114 : }
371 :
372 : //----------------------------------------------------------------------
373 : // Odds are that it changed.....
374 338 : Bool SDGrid::changed(const VisBuffer& /*vb*/) {
375 338 : return false;
376 : }
377 :
378 : //----------------------------------------------------------------------
379 0 : SDGrid::SDGrid(const SDGrid& other)
380 : : FTMachine(),
381 0 : rowPixel {MaskedPixelRef(xyPos)},
382 0 : cache {Cache(*const_cast<SDGrid *>(this))}
383 : {
384 0 : operator=(other);
385 0 : }
386 :
387 : #define NEED_UNDERSCORES
388 : #if defined(NEED_UNDERSCORES)
389 : #define grdsf grdsf_
390 : #define grdgauss grdgauss_
391 : #define grdjinc1 grdjinc1_
392 : #endif
393 :
394 : extern "C" {
395 : void grdsf(Double*, Double*);
396 : void grdgauss(Double*, Double*, Double*);
397 : void grdjinc1(Double*, Double*, Int*, Double*);
398 : }
399 :
400 0 : void SDGrid::dumpConvolutionFunction(const casacore::String &outfile,
401 : const casacore::Vector<casacore::Float> &f) const {
402 :
403 0 : ofstream ofs(outfile.c_str());
404 0 : for (size_t i = 0 ; i < f.nelements() ; i++) {
405 0 : ofs << i << " " << f[i] << endl;
406 : }
407 0 : ofs.close();
408 0 : }
409 :
410 : //----------------------------------------------------------------------
411 238 : void SDGrid::init() {
412 :
413 238 : logIO() << LogOrigin("SDGrid", "init") << LogIO::NORMAL;
414 :
415 238 : ok();
416 :
417 238 : isTiled = false;
418 238 : nx = image->shape()(0);
419 238 : ny = image->shape()(1);
420 238 : npol = image->shape()(2);
421 238 : nchan = image->shape()(3);
422 :
423 238 : sumWeight.resize(npol, nchan);
424 :
425 : // Set up image cache needed for gridding. For BOX-car convolution
426 : // we can use non-overlapped tiles. Otherwise we need to use
427 : // overlapped tiles and additive gridding so that only increments
428 : // to a tile are written.
429 238 : if (imageCache) delete imageCache;
430 238 : imageCache = nullptr;
431 :
432 238 : convType = downcase(convType);
433 238 : logIO() << "Convolution function : " << convType << LogIO::DEBUG1 << LogIO::POST;
434 :
435 : // Compute convolution function
436 238 : if (convType=="pb") {
437 : // Do nothing
438 : }
439 192 : else if (convType=="box") {
440 174 : convSupport=(userSetSupport_p >= 0) ? userSetSupport_p : 0;
441 174 : logIO() << "Support : " << convSupport << " pixels" << LogIO::POST;
442 174 : convSampling=100;
443 174 : convSize=convSampling*(2*convSupport+2);
444 174 : convFunc.resize(convSize);
445 174 : convFunc=0.0;
446 17574 : for (Int i=0;i<convSize/2;i++) {
447 17400 : convFunc(i)=1.0;
448 : }
449 : }
450 18 : else if (convType=="sf") {
451 12 : convSupport=(userSetSupport_p >= 0) ? userSetSupport_p : 3;
452 12 : logIO() << "Support : " << convSupport << " pixels" << LogIO::POST;
453 12 : convSampling=100;
454 12 : convSize=convSampling*(2*convSupport+2);
455 12 : convFunc.resize(convSize);
456 12 : convFunc=0.0;
457 6612 : for (Int i=0;i<convSampling*convSupport;i++) {
458 6600 : Double nu=Double(i)/Double(convSupport*convSampling);
459 : Double val;
460 6600 : grdsf(&nu, &val);
461 6600 : convFunc(i)=(1.0-nu*nu)*val;
462 : }
463 : }
464 6 : else if (convType=="gauss") {
465 : // default is b=1.0 (Mangum et al. 2007)
466 2 : Double hwhm=(gwidth_p > 0.0) ? Double(gwidth_p) : sqrt(log(2.0));
467 2 : Float truncate=(truncate_p >= 0.0) ? truncate_p : 3.0*hwhm;
468 2 : convSampling=100;
469 2 : Int itruncate=(Int)(truncate*Double(convSampling)+0.5);
470 2 : logIO() << LogIO::DEBUG1 << "hwhm=" << hwhm << LogIO::POST;
471 2 : logIO() << LogIO::DEBUG1 << "itruncate=" << itruncate << LogIO::POST;
472 : //convSupport=(Int)(truncate+0.5);
473 2 : convSupport = (Int)(truncate);
474 2 : convSupport += (((truncate-(Float)convSupport) > 0.0) ? 1 : 0);
475 2 : convSize=convSampling*(2*convSupport+2);
476 2 : convFunc.resize(convSize);
477 2 : convFunc=0.0;
478 : Double val, x;
479 504 : for (Int i = 0 ; i <= itruncate ; i++) {
480 502 : x = Double(i)/Double(convSampling);
481 502 : grdgauss(&hwhm, &x, &val);
482 502 : convFunc(i)=val;
483 : }
484 : }
485 4 : else if (convType=="gjinc") {
486 : // default is b=2.52, c=1.55 (Mangum et al. 2007)
487 4 : Double hwhm = (gwidth_p > 0.0) ? Double(gwidth_p) : sqrt(log(2.0))*2.52;
488 4 : Double c = (jwidth_p > 0.0) ? Double(jwidth_p) : 1.55;
489 : //Float truncate = truncate_p;
490 4 : convSampling = 100;
491 4 : Int itruncate=(Int)(truncate_p*Double(convSampling)+0.5);
492 4 : logIO() << LogIO::DEBUG1 << "hwhm=" << hwhm << LogIO::POST;
493 4 : logIO() << LogIO::DEBUG1 << "c=" << c << LogIO::POST;
494 4 : logIO() << LogIO::DEBUG1 << "itruncate=" << itruncate << LogIO::POST;
495 : //convSupport=(truncate_p >= 0.0) ? (Int)(truncate_p+0.5) : (Int)(2*c+0.5);
496 4 : Float convSupportF = (truncate_p >= 0.0) ? truncate_p : (2*c);
497 4 : convSupport = (Int)convSupportF;
498 4 : convSupport += (((convSupportF-(Float)convSupport) > 0.0) ? 1 : 0);
499 4 : convSize=convSampling*(2*convSupport+2);
500 4 : convFunc.resize(convSize);
501 4 : convFunc=0.0;
502 : //UNUSED: Double r;
503 : Double x, val1, val2;
504 4 : Int normalize = 1;
505 4 : if (itruncate >= 0) {
506 0 : for (Int i = 0 ; i < itruncate ; i++) {
507 0 : x = Double(i) / Double(convSampling);
508 : //r = Double(i) / (Double(hwhm)*Double(convSampling));
509 0 : grdgauss(&hwhm, &x, &val1);
510 0 : grdjinc1(&c, &x, &normalize, &val2);
511 0 : convFunc(i) = val1 * val2;
512 : }
513 : }
514 : else {
515 : // default is to truncate at first null
516 764 : for (Int i=0;i<convSize;i++) {
517 764 : x = Double(i) / Double(convSampling);
518 : //r = Double(i) / (Double(hwhm)*Double(convSampling));
519 764 : grdjinc1(&c, &x, &normalize, &val2);
520 764 : if (val2 <= 0.0) {
521 : logIO() << LogIO::DEBUG1
522 : << "convFunc is automatically truncated at radius "
523 4 : << x << LogIO::POST;
524 4 : break;
525 : }
526 760 : grdgauss(&hwhm, &x, &val1);
527 760 : convFunc(i) = val1 * val2;
528 : }
529 : }
530 : }
531 : else {
532 0 : logIO_p << "Unknown convolution function " << convType << LogIO::EXCEPTION;
533 : }
534 :
535 : // Convolution function debug
536 : // String outfile = convType + ".dat";
537 : // dumpConvolutionFunction(outfile,convFunc);
538 :
539 238 : if (wImage) delete wImage;
540 238 : wImage = new TempImage<Float>(image->shape(), image->coordinates());
541 :
542 238 : }
543 :
544 124 : void SDGrid::collectPerfs(){
545 : #if defined(SDGRID_PERFS)
546 372 : LogIO os(LogOrigin(name(), "collectPerfs"));
547 : std::vector<std::string> probes_json {
548 : cNextChunk.json()
549 : , cMatchAllSpwChans.json()
550 : , cMatchChannel.json()
551 : , cPickWeights.json()
552 : , cInterpolateFrequencyToGrid.json()
553 : , cSearchValidPointing.json()
554 : , cComputeSplines.json()
555 : , cResetFrame.json()
556 : , cInterpolateDirection.json()
557 : , cConvertDirection.json()
558 : , cComputeDirectionPixel.json()
559 : , cHandleMovingSource.json()
560 : , cGridData.json()
561 1860 : };
562 : os << "PERFS<SDGRID> "
563 : << "{ \"note\": \"sum, min, mean, max are in units of nanoseconds.\""
564 : << ", \"probes\": "
565 248 : << "{ " << join(probes_json.data(), probes_json.size(), ", " ) << " }"
566 : << " }"
567 248 : << LogIO::POST;
568 : #endif
569 124 : }
570 :
571 :
572 : // This is nasty, we should use CountedPointers here.
573 248 : SDGrid::~SDGrid() {
574 : //fclose(pfile);
575 124 : if (imageCache) delete imageCache; imageCache = 0;
576 124 : if (arrayLattice) delete arrayLattice; arrayLattice = 0;
577 124 : if (wImage) delete wImage; wImage = 0;
578 124 : if (wImageCache) delete wImageCache; wImageCache = 0;
579 124 : if (wArrayLattice) delete wArrayLattice; wArrayLattice = 0;
580 124 : if (interpolator) delete interpolator; interpolator = 0;
581 :
582 124 : collectPerfs();
583 248 : }
584 :
585 46 : void SDGrid::findPBAsConvFunction(const ImageInterface<Complex>& image,
586 : const VisBuffer& vb) {
587 :
588 : // Get the coordinate system and increase the sampling by
589 : // a factor of ~ 100.
590 92 : CoordinateSystem coords(image.coordinates());
591 :
592 : // Set up the convolution function: make the buffer plenty
593 : // big so that we can trim it back
594 46 : convSupport=max(128, sj_p->support(vb, coords));
595 :
596 46 : convSampling=100;
597 46 : convSize=convSampling*convSupport;
598 :
599 : // Make a one dimensional image to calculate the
600 : // primary beam. We oversample this by a factor of
601 : // convSampling.
602 46 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
603 46 : AlwaysAssert(directionIndex>=0, AipsError);
604 92 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
605 92 : Vector<Double> sampling;
606 46 : sampling = dc.increment();
607 46 : sampling*=1.0/Double(convSampling);
608 46 : dc.setIncrement(sampling);
609 :
610 : // Set the reference value to the first pointing in the coordinate
611 : // system used in the POINTING table.
612 : {
613 46 : uInt row=0;
614 :
615 : // reset lastAntID_p to use correct antenna position
616 46 : lastAntID_p = -1;
617 :
618 46 : const MSPointingColumns& act_mspc = vb.msColumns().pointing();
619 46 : Bool nullPointingTable=(act_mspc.nrow() < 1);
620 : // uInt pointIndex=getIndex(*mspc, vb.time()(row), vb.timeInterval()(row), vb.antenna1()(row));
621 46 : Int pointIndex=-1;
622 46 : if(!nullPointingTable){
623 : //if(vb.newMS()) This thing is buggy...using msId change
624 46 : if(vb.msId() != msId_p){
625 28 : lastIndex_p=0;
626 28 : if(lastIndexPerAnt_p.nelements() < (size_t)vb.numberAnt()) {
627 28 : lastIndexPerAnt_p.resize(vb.numberAnt());
628 : }
629 28 : lastIndexPerAnt_p=0;
630 28 : msId_p=vb.msId();
631 : }
632 46 : pointIndex=getIndex(act_mspc, vb.time()(row), -1.0, vb.antenna1()(row));
633 46 : if(pointIndex < 0)
634 0 : pointIndex=getIndex(act_mspc, vb.time()(row), vb.timeInterval()(row), vb.antenna1()(row));
635 :
636 : }
637 46 : if(!nullPointingTable && ((pointIndex<0)||(pointIndex>=Int(act_mspc.time().nrow())))) {
638 0 : ostringstream o;
639 0 : o << "Failed to find pointing information for time " <<
640 0 : MVTime(vb.time()(row)/86400.0);
641 0 : logIO_p << String(o) << LogIO::EXCEPTION;
642 : }
643 92 : MEpoch epoch(Quantity(vb.time()(row), "s"));
644 46 : if(!pointingToImage) {
645 46 : lastAntID_p=vb.antenna1()(row);
646 92 : MPosition pos;
647 46 : pos=vb.msColumns().antenna().positionMeas()(lastAntID_p);
648 46 : mFrame_p=MeasFrame(epoch, pos);
649 46 : if(!nullPointingTable){
650 46 : worldPosMeas=directionMeas(act_mspc, pointIndex);
651 : }
652 : else{
653 0 : worldPosMeas=vb.direction1()(row);
654 : }
655 : // Make a machine to convert from the worldPosMeas to the output
656 : // Direction Measure type for the relevant frame
657 92 : MDirection::Ref outRef(dc.directionType(), mFrame_p);
658 46 : pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
659 :
660 46 : if(!pointingToImage) {
661 0 : logIO_p << "Cannot make direction conversion machine" << LogIO::EXCEPTION;
662 : }
663 : }
664 : else {
665 0 : mFrame_p.resetEpoch(epoch);
666 0 : if(lastAntID_p != vb.antenna1()(row)){
667 0 : logIO_p << LogIO::DEBUGGING
668 : << "updating antenna position. MS ID " << msId_p
669 : << ", last antenna ID " << lastAntID_p
670 0 : << " new antenna ID " << vb.antenna1()(row) << LogIO::POST;
671 0 : MPosition pos;
672 0 : lastAntID_p=vb.antenna1()(row);
673 0 : pos=vb.msColumns().antenna().positionMeas()(lastAntID_p);
674 0 : mFrame_p.resetPosition(pos);
675 : }
676 : }
677 46 : if(!nullPointingTable){
678 46 : worldPosMeas=(*pointingToImage)(directionMeas(act_mspc,pointIndex));
679 : }
680 : else{
681 0 : worldPosMeas=(*pointingToImage)(vb.direction1()(row));
682 : }
683 46 : delete pointingToImage;
684 46 : pointingToImage=0;
685 : }
686 :
687 46 : directionCoord=coords.directionCoordinate(directionIndex);
688 : //make sure we use the same units
689 46 : worldPosMeas.set(dc.worldAxisUnits()(0));
690 :
691 : // Reference pixel may be modified in dc.setReferenceValue when
692 : // projection type is SFL. To take into account this effect,
693 : // keep original reference pixel here and subtract it from
694 : // the reference pixel after dc.setReferenceValue instead
695 : // of setting reference pixel to (0,0).
696 92 : Vector<Double> const originalReferencePixel = dc.referencePixel();
697 46 : dc.setReferenceValue(worldPosMeas.getAngle().getValue());
698 : //Vector<Double> unitVec(2);
699 : //unitVec=0.0;
700 : //dc.setReferencePixel(unitVec);
701 138 : Vector<Double> updatedReferencePixel = dc.referencePixel() - originalReferencePixel;
702 46 : dc.setReferencePixel(updatedReferencePixel);
703 :
704 46 : coords.replaceCoordinate(dc, directionIndex);
705 :
706 92 : IPosition pbShape(4, convSize, 2, 1, 1);
707 92 : IPosition start(4, 0, 0, 0, 0);
708 :
709 92 : TempImage<Complex> onedPB(pbShape, coords);
710 :
711 46 : onedPB.set(Complex(1.0, 0.0));
712 :
713 46 : AlwaysAssert(sj_p, AipsError);
714 46 : sj_p->apply(onedPB, onedPB, vb, 0);
715 :
716 92 : IPosition pbSlice(4, convSize, 1, 1, 1);
717 138 : Vector<Float> tempConvFunc=real(onedPB.getSlice(start, pbSlice, true));
718 : // Find number of significant points
719 46 : uInt cfLen=0;
720 168404 : for(uInt i=0;i<tempConvFunc.nelements();++i) {
721 168404 : if(tempConvFunc(i)<1e-3) break;
722 168358 : ++cfLen;
723 : }
724 46 : if(cfLen<1) {
725 : logIO() << LogIO::SEVERE
726 : << "Possible problem in primary beam calculation: no points in gridding function"
727 0 : << " - no points to be gridded on this image?" << LogIO::POST;
728 0 : cfLen=1;
729 : }
730 46 : Vector<Float> trimConvFunc=tempConvFunc(Slice(0,cfLen-1,1));
731 :
732 : // Now fill in the convolution function vector
733 46 : convSupport=cfLen/convSampling;
734 46 : convSize=convSampling*(2*convSupport+2);
735 46 : convFunc.resize(convSize);
736 46 : convFunc=0.0;
737 46 : convFunc(Slice(0,cfLen-1,1))=trimConvFunc(Slice(0,cfLen-1,1));
738 :
739 :
740 46 : }
741 :
742 : // Initialize for a transform from the Sky domain. This means that
743 : // we grid-correct, and FFT the image
744 10 : void SDGrid::initializeToVis(ImageInterface<Complex>& iimage,
745 : const VisBuffer& vb)
746 : {
747 10 : image=&iimage;
748 :
749 10 : ok();
750 :
751 10 : init();
752 :
753 10 : if(convType=="pb") {
754 10 : findPBAsConvFunction(*image, vb);
755 : }
756 :
757 : // reset msId_p and lastAntID_p to -1
758 : // this is to ensure correct antenna position in getXYPos
759 10 : msId_p = -1;
760 10 : lastAntID_p = -1;
761 :
762 : // Initialize the maps for polarization and channel. These maps
763 : // translate visibility indices into image indices
764 10 : initMaps(vb);
765 :
766 : // First get the CoordinateSystem for the image and then find
767 : // the DirectionCoordinate
768 20 : CoordinateSystem coords=image->coordinates();
769 10 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
770 10 : AlwaysAssert(directionIndex>=0, AipsError);
771 10 : directionCoord=coords.directionCoordinate(directionIndex);
772 : /*if((image->shape().product())>cachesize) {
773 : isTiled=true;
774 : }
775 : else {
776 : isTiled=false;
777 : }*/
778 10 : isTiled=false;
779 10 : nx = image->shape()(0);
780 10 : ny = image->shape()(1);
781 10 : npol = image->shape()(2);
782 10 : nchan = image->shape()(3);
783 :
784 : // If we are memory-based then read the image in and create an
785 : // ArrayLattice otherwise just use the PagedImage
786 : /*if(isTiled) {
787 : lattice=image;
788 : wLattice=wImage;
789 : }
790 : else*/
791 : {
792 : // Make the grid the correct shape and turn it into an array lattice
793 20 : IPosition gridShape(4, nx, ny, npol, nchan);
794 10 : griddedData.resize(gridShape);
795 10 : griddedData = Complex(0.0);
796 :
797 10 : wGriddedData.resize(gridShape);
798 10 : wGriddedData = 0.0;
799 :
800 10 : if(arrayLattice) delete arrayLattice; arrayLattice=0;
801 10 : arrayLattice = new ArrayLattice<Complex>(griddedData);
802 :
803 10 : if(wArrayLattice) delete wArrayLattice; wArrayLattice=0;
804 10 : wArrayLattice = new ArrayLattice<Float>(wGriddedData);
805 10 : wArrayLattice->set(0.0);
806 10 : wLattice=wArrayLattice;
807 :
808 : // Now find the SubLattice corresponding to the image
809 30 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
810 20 : IPosition stride(4, 1);
811 30 : IPosition trc(blc+image->shape()-stride);
812 20 : LCBox gridBox(blc, trc, gridShape);
813 20 : SubLattice<Complex> gridSub(*arrayLattice, gridBox, true);
814 :
815 : // Do the copy
816 10 : gridSub.copyData(*image);
817 :
818 10 : lattice=arrayLattice;
819 : }
820 :
821 10 : AlwaysAssert(lattice, AipsError);
822 10 : AlwaysAssert(wLattice, AipsError);
823 :
824 10 : }
825 :
826 0 : void SDGrid::finalizeToVis()
827 : {
828 : /*if(isTiled) {
829 :
830 : logIO() << LogOrigin("SDGrid", "finalizeToVis") << LogIO::NORMAL;
831 :
832 : AlwaysAssert(imageCache, AipsError);
833 : AlwaysAssert(image, AipsError);
834 : ostringstream o;
835 : imageCache->flush();
836 : imageCache->showCacheStatistics(o);
837 : logIO() << o.str() << LogIO::POST;
838 : }*/
839 0 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
840 0 : }
841 :
842 :
843 : // Initialize the FFT to the Sky.
844 : // Here we have to setup and initialize the grid.
845 228 : void SDGrid::initializeToSky(ImageInterface<Complex>& iimage,
846 : Matrix<Float>& weight, const VisBuffer& vb)
847 : {
848 : // image always points to the image
849 228 : image=&iimage;
850 :
851 228 : ok();
852 :
853 228 : init();
854 :
855 228 : if (convType == "pb") findPBAsConvFunction(*image, vb);
856 :
857 : // reset msId_p and lastAntID_p to -1
858 : // this is to ensure correct antenna position in getXYPos
859 228 : msId_p = -1;
860 228 : lastAntID_p = -1;
861 :
862 : // Initialize the maps for polarization and channel.
863 : // These maps translate visibility indices into image indices
864 228 : initMaps(vb);
865 :
866 : // No longer using isTiled
867 228 : isTiled=false;
868 228 : nx = image->shape()(0);
869 228 : ny = image->shape()(1);
870 228 : npol = image->shape()(2);
871 228 : nchan = image->shape()(3);
872 :
873 228 : sumWeight = 0.0;
874 228 : weight.resize(sumWeight.shape());
875 228 : weight = 0.0;
876 :
877 : // First get the CoordinateSystem for the image
878 : // and then find the DirectionCoordinate
879 456 : CoordinateSystem coords = image->coordinates();
880 228 : Int directionIndex = coords.findCoordinate(Coordinate::DIRECTION);
881 228 : AlwaysAssert(directionIndex >= 0, AipsError);
882 228 : directionCoord = coords.directionCoordinate(directionIndex);
883 :
884 : // Initialize for in memory gridding.
885 : // lattice will point to the ArrayLattice
886 456 : IPosition gridShape(4, nx, ny, npol, nchan);
887 456 : logIO() << LogOrigin("SDGrid", "initializeToSky", WHERE) << LogIO::DEBUGGING
888 456 : << "gridShape = " << gridShape << LogIO::POST;
889 :
890 228 : griddedData.resize(gridShape);
891 228 : griddedData = Complex(0.0);
892 228 : if (arrayLattice) delete arrayLattice;
893 228 : arrayLattice = new ArrayLattice<Complex>(griddedData);
894 228 : lattice = arrayLattice;
895 228 : AlwaysAssert(lattice, AipsError);
896 :
897 228 : wGriddedData.resize(gridShape);
898 228 : wGriddedData=0.0;
899 228 : if (wArrayLattice) delete wArrayLattice;
900 228 : wArrayLattice = new ArrayLattice<Float>(wGriddedData);
901 228 : wLattice = wArrayLattice;
902 228 : AlwaysAssert(wLattice, AipsError);
903 :
904 : // clipping related stuff
905 228 : if (clipminmax_) {
906 16 : gmin_.resize(gridShape);
907 16 : gmin_ = Complex(FLT_MAX);
908 16 : gmax_.resize(gridShape);
909 16 : gmax_ = Complex(-FLT_MAX);
910 16 : wmin_.resize(gridShape);
911 16 : wmin_ = 0.0f;
912 16 : wmax_.resize(gridShape);
913 16 : wmax_ = 0.0f;
914 16 : npoints_.resize(gridShape);
915 16 : npoints_ = 0;
916 : }
917 :
918 : // debug messages
919 684 : LogOrigin msgOrigin("SDGrid", "initializeToSky", WHERE);
920 228 : auto & logger = logIO();
921 228 : logger << msgOrigin << LogIO::DEBUGGING;
922 228 : logger.output() << "clipminmax_ = " << std::boolalpha << clipminmax_
923 228 : << " (" << std::noboolalpha << clipminmax_ << ")";
924 228 : logger << LogIO::POST;
925 228 : if (clipminmax_) {
926 : logger << msgOrigin << LogIO::DEBUGGING
927 : << "will use clipping-capable Fortran gridder ggridsd2 for imaging"
928 16 : << LogIO::POST;
929 : }
930 228 : }
931 :
932 228 : void SDGrid::finalizeToSky()
933 : {
934 228 : if (pointingToImage) delete pointingToImage;
935 228 : pointingToImage = nullptr;
936 228 : }
937 :
938 0 : Array<Complex>* SDGrid::getDataPointer(const IPosition& centerLoc2D,
939 : Bool readonly) {
940 : Array<Complex>* result;
941 : // Is tiled: get tiles and set up offsets
942 0 : centerLoc(0)=centerLoc2D(0);
943 0 : centerLoc(1)=centerLoc2D(1);
944 0 : result=&imageCache->tile(offsetLoc, centerLoc, readonly);
945 0 : return result;
946 : }
947 0 : Array<Float>* SDGrid::getWDataPointer(const IPosition& centerLoc2D,
948 : Bool readonly) {
949 : Array<Float>* result;
950 : // Is tiled: get tiles and set up offsets
951 0 : centerLoc(0)=centerLoc2D(0);
952 0 : centerLoc(1)=centerLoc2D(1);
953 0 : result=&wImageCache->tile(offsetLoc, centerLoc, readonly);
954 0 : return result;
955 : }
956 :
957 : #define NEED_UNDERSCORES
958 : #if defined(NEED_UNDERSCORES)
959 : #define ggridsd ggridsd_
960 : #define dgridsd dgridsd_
961 : #define ggridsdclip ggridsdclip_
962 : #endif
963 :
964 : extern "C" {
965 : void ggridsd(Double*,
966 : const Complex*,
967 : Int*,
968 : Int*,
969 : Int*,
970 : const Int*,
971 : const Int*,
972 : const Float*,
973 : Int*,
974 : Int*,
975 : Complex*,
976 : Float*,
977 : Int*,
978 : Int*,
979 : Int *,
980 : Int *,
981 : Int*,
982 : Int*,
983 : Float*,
984 : Int*,
985 : Int*,
986 : Double*);
987 : void ggridsdclip(Double*,
988 : const Complex*,
989 : Int*,
990 : Int*,
991 : Int*,
992 : const Int*,
993 : const Int*,
994 : const Float*,
995 : Int*,
996 : Int*,
997 : Complex*,
998 : Float*,
999 : Int*,
1000 : Complex*,
1001 : Float*,
1002 : Complex*,
1003 : Float*,
1004 : Int*,
1005 : Int*,
1006 : Int *,
1007 : Int *,
1008 : Int*,
1009 : Int*,
1010 : Float*,
1011 : Int*,
1012 : Int*,
1013 : Double*);
1014 : void dgridsd(Double*,
1015 : Complex*,
1016 : Int*,
1017 : Int*,
1018 : const Int*,
1019 : const Int*,
1020 : Int*,
1021 : Int*,
1022 : const Complex*,
1023 : Int*,
1024 : Int*,
1025 : Int *,
1026 : Int *,
1027 : Int*,
1028 : Int*,
1029 : Float*,
1030 : Int*,
1031 : Int*);
1032 : }
1033 :
1034 4322 : void SDGrid::put(const VisBuffer& vb, Int row, Bool dopsf,
1035 : FTMachine::Type type)
1036 : {
1037 8644 : LogIO os(LogOrigin("SDGrid", "put"));
1038 :
1039 4322 : gridOk(convSupport);
1040 :
1041 : // Check if ms has changed then cache new spw and chan selection
1042 4322 : if (vb.newMS()) {
1043 : #if defined(SDGRID_PERFS)
1044 43 : StartStop trigger(cMatchAllSpwChans);
1045 : #endif
1046 43 : matchAllSpwChans(vb);
1047 43 : lastIndex_p = 0;
1048 43 : if (lastIndexPerAnt_p.nelements() < (size_t)vb.numberAnt()) {
1049 0 : lastIndexPerAnt_p.resize(vb.numberAnt());
1050 : }
1051 43 : lastIndexPerAnt_p = 0;
1052 : }
1053 :
1054 : // Here we redo the match or use previous match
1055 : // Channel matching for the actual spectral window of buffer
1056 4322 : if (doConversion_p[vb.spectralWindow()]) {
1057 : #if defined(SDGRID_PERFS)
1058 964 : StartStop trigger(cMatchChannel);
1059 : #endif
1060 964 : matchChannel(vb.spectralWindow(), vb);
1061 : }
1062 : else {
1063 3358 : chanMap.resize();
1064 3358 : chanMap = multiChanMap_p[vb.spectralWindow()];
1065 : }
1066 :
1067 : // No point in reading data if its not matching in frequency
1068 4322 : if (max(chanMap)==-1)
1069 0 : return;
1070 :
1071 8644 : Matrix<Float> imagingweight;
1072 4322 : pickWeights(vb, imagingweight);
1073 :
1074 4322 : if (type==FTMachine::PSF || type==FTMachine::COVERAGE)
1075 2161 : dopsf = true;
1076 4322 : if (dopsf) type=FTMachine::PSF;
1077 :
1078 8644 : Cube<Complex> data;
1079 : //Fortran gridder need the flag as ints
1080 8644 : Cube<Int> flags;
1081 8644 : Matrix<Float> elWeight;
1082 : {
1083 : #if defined(SDGRID_PERFS)
1084 8644 : StartStop trigger(cInterpolateFrequencyToGrid);
1085 : #endif
1086 4322 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
1087 : }
1088 :
1089 : Bool iswgtCopy;
1090 : const Float *wgtStorage;
1091 4322 : wgtStorage=elWeight.getStorage(iswgtCopy);
1092 : Bool isCopy;
1093 4322 : const Complex *datStorage = nullptr;
1094 4322 : if (!dopsf)
1095 2161 : datStorage = data.getStorage(isCopy);
1096 :
1097 : // If row is -1 then we pass through all rows
1098 : Int startRow, endRow, nRow;
1099 4322 : if (row == -1) {
1100 4322 : nRow = vb.nRow();
1101 4322 : startRow = 0;
1102 4322 : endRow = nRow - 1;
1103 : } else {
1104 0 : nRow = 1;
1105 0 : startRow = endRow = row;
1106 : }
1107 :
1108 8644 : Vector<Int> rowFlags(vb.flagRow().nelements(),0);
1109 356400 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1110 352078 : if (vb.flagRow()(rownr)) rowFlags(rownr) = 1;
1111 : }
1112 :
1113 : // Take care of translation of Bools to Integer
1114 4322 : Int idopsf = dopsf ? 1 : 0;
1115 :
1116 : /*if(isTiled) {
1117 : }
1118 : else*/
1119 : {
1120 8644 : Matrix<Double> xyPositions(2, endRow - startRow + 1);
1121 4322 : xyPositions = -1e9; // make sure failed getXYPos does not fall on grid
1122 356400 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1123 352078 : if (getXYPos(vb, rownr)) {
1124 352078 : xyPositions(0, rownr) = xyPos(0);
1125 352078 : xyPositions(1, rownr) = xyPos(1);
1126 : }
1127 : }
1128 : {
1129 : #if defined(SDGRID_PERFS)
1130 8644 : StartStop trigger(cGridData);
1131 : #endif
1132 : Bool del;
1133 : // IPosition s(data.shape());
1134 4322 : const IPosition& fs=flags.shape();
1135 8644 : std::vector<Int> s(fs.begin(), fs.end());
1136 : Bool datCopy, wgtCopy;
1137 4322 : Complex * datStor=griddedData.getStorage(datCopy);
1138 4322 : Float * wgtStor=wGriddedData.getStorage(wgtCopy);
1139 :
1140 4322 : Bool call_ggridsd = !clipminmax_ || dopsf;
1141 :
1142 4322 : if (call_ggridsd) {
1143 :
1144 8620 : ggridsd(xyPositions.getStorage(del),
1145 : datStorage,
1146 4310 : &s[0],
1147 4310 : &s[1],
1148 : &idopsf,
1149 4310 : flags.getStorage(del),
1150 4310 : rowFlags.getStorage(del),
1151 : wgtStorage,
1152 4310 : &s[2],
1153 : &row,
1154 : datStor,
1155 : wgtStor,
1156 : &nx,
1157 : &ny,
1158 : &npol,
1159 : &nchan,
1160 : &convSupport,
1161 : &convSampling,
1162 : convFunc.getStorage(del),
1163 : chanMap.getStorage(del),
1164 : polMap.getStorage(del),
1165 : sumWeight.getStorage(del));
1166 :
1167 : }
1168 : else {
1169 : Bool gminCopy;
1170 12 : Complex *gminStor = gmin_.getStorage(gminCopy);
1171 : Bool gmaxCopy;
1172 12 : Complex *gmaxStor = gmax_.getStorage(gmaxCopy);
1173 : Bool wminCopy;
1174 12 : Float *wminStor = wmin_.getStorage(wminCopy);
1175 : Bool wmaxCopy;
1176 12 : Float *wmaxStor = wmax_.getStorage(wmaxCopy);
1177 : Bool npCopy;
1178 12 : Int *npStor = npoints_.getStorage(npCopy);
1179 :
1180 24 : ggridsdclip(xyPositions.getStorage(del),
1181 : datStorage,
1182 12 : &s[0],
1183 12 : &s[1],
1184 : &idopsf,
1185 12 : flags.getStorage(del),
1186 12 : rowFlags.getStorage(del),
1187 : wgtStorage,
1188 12 : &s[2],
1189 : &row,
1190 : datStor,
1191 : wgtStor,
1192 : npStor,
1193 : gminStor,
1194 : wminStor,
1195 : gmaxStor,
1196 : wmaxStor,
1197 : &nx,
1198 : &ny,
1199 : &npol,
1200 : &nchan,
1201 : &convSupport,
1202 : &convSampling,
1203 : convFunc.getStorage(del),
1204 : chanMap.getStorage(del),
1205 : polMap.getStorage(del),
1206 : sumWeight.getStorage(del));
1207 :
1208 12 : gmin_.putStorage(gminStor, gminCopy);
1209 12 : gmax_.putStorage(gmaxStor, gmaxCopy);
1210 12 : wmin_.putStorage(wminStor, wminCopy);
1211 12 : wmax_.putStorage(wmaxStor, wmaxCopy);
1212 12 : npoints_.putStorage(npStor, npCopy);
1213 : }
1214 4322 : griddedData.putStorage(datStor, datCopy);
1215 4322 : wGriddedData.putStorage(wgtStor, wgtCopy);
1216 : }
1217 : }
1218 4322 : if (!dopsf)
1219 2161 : data.freeStorage(datStorage, isCopy);
1220 4322 : elWeight.freeStorage(wgtStorage,iswgtCopy);
1221 : }
1222 :
1223 338 : void SDGrid::get(VisBuffer& vb, Int row)
1224 : {
1225 676 : LogIO os(LogOrigin("SDGrid", "get"));
1226 :
1227 338 : gridOk(convSupport);
1228 :
1229 : // If row is -1 then we pass through all rows
1230 : Int startRow, endRow, nRow;
1231 338 : if (row==-1) {
1232 338 : nRow=vb.nRow();
1233 338 : startRow=0;
1234 338 : endRow=nRow-1;
1235 338 : vb.modelVisCube()=Complex(0.0,0.0);
1236 : } else {
1237 0 : nRow=1;
1238 0 : startRow=row;
1239 0 : endRow=row;
1240 0 : vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
1241 : }
1242 :
1243 :
1244 : //Check if ms has changed then cache new spw and chan selection
1245 338 : if(vb.newMS()){
1246 0 : matchAllSpwChans(vb);
1247 0 : lastIndex_p=0;
1248 0 : if (lastIndexPerAnt_p.nelements() < (size_t)vb.numberAnt()) {
1249 0 : lastIndexPerAnt_p.resize(vb.numberAnt());
1250 : }
1251 0 : lastIndexPerAnt_p=0;
1252 : }
1253 :
1254 : //Here we redo the match or use previous match
1255 :
1256 : //Channel matching for the actual spectral window of buffer
1257 338 : if(doConversion_p[vb.spectralWindow()]){
1258 0 : matchChannel(vb.spectralWindow(), vb);
1259 : }
1260 : else{
1261 338 : chanMap.resize();
1262 338 : chanMap=multiChanMap_p[vb.spectralWindow()];
1263 : }
1264 :
1265 : //No point in reading data if its not matching in frequency
1266 338 : if(max(chanMap)==-1)
1267 0 : return;
1268 676 : Cube<Complex> data;
1269 676 : Cube<Int> flags;
1270 338 : getInterpolateArrays(vb, data, flags);
1271 :
1272 : Complex *datStorage;
1273 : Bool isCopy;
1274 338 : datStorage=data.getStorage(isCopy);
1275 : // NOTE: with MS V2.0 the pointing could change per antenna and timeslot
1276 : //
1277 676 : Vector<Int> rowFlags(vb.flagRow().nelements());
1278 338 : rowFlags=0;
1279 3346 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1280 3008 : if(vb.flagRow()(rownr)) rowFlags(rownr)=1;
1281 : //single dish yes ?
1282 3008 : if(max(vb.uvw()(rownr).vector()) > 0.0) rowFlags(rownr)=1;
1283 : }
1284 :
1285 :
1286 : /*if(isTiled) {
1287 :
1288 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1289 :
1290 : if(getXYPos(vb, rownr)) {
1291 :
1292 : // Get the tile
1293 : IPosition centerLoc2D(2, Int(xyPos(0)), Int(xyPos(1)));
1294 : Array<Complex>* dataPtr=getDataPointer(centerLoc2D, true);
1295 : Int aNx=dataPtr->shape()(0);
1296 : Int aNy=dataPtr->shape()(1);
1297 :
1298 : // Now use FORTRAN to do the gridding. Remember to
1299 : // ensure that the shape and offsets of the tile are
1300 : // accounted for.
1301 : Bool del;
1302 : Vector<Double> actualPos(2);
1303 : for (Int i=0;i<2;i++) {
1304 : actualPos(i)=xyPos(i)-Double(offsetLoc(i));
1305 : }
1306 : // IPosition s(data.shape());
1307 : const IPosition& fs=data.shape();
1308 : std::vector<Int> s(fs.begin(), fs.end());
1309 :
1310 : dgridsd(actualPos.getStorage(del),
1311 : datStorage,
1312 : &s[0],
1313 : &s[1],
1314 : flags.getStorage(del),
1315 : rowFlags.getStorage(del),
1316 : &s[2],
1317 : &rownr,
1318 : dataPtr->getStorage(del),
1319 : &aNx,
1320 : &aNy,
1321 : &npol,
1322 : &nchan,
1323 : &convSupport,
1324 : &convSampling,
1325 : convFunc.getStorage(del),
1326 : chanMap.getStorage(del),
1327 : polMap.getStorage(del));
1328 : }
1329 : }
1330 : }
1331 : else*/
1332 : {
1333 676 : Matrix<Double> xyPositions(2, endRow-startRow+1);
1334 338 : xyPositions=-1e9;
1335 3346 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1336 3008 : if(getXYPos(vb, rownr)) {
1337 3008 : xyPositions(0, rownr)=xyPos(0);
1338 3008 : xyPositions(1, rownr)=xyPos(1);
1339 : }
1340 : }
1341 :
1342 : Bool del;
1343 : // IPosition s(data.shape());
1344 338 : const IPosition& fs=data.shape();
1345 676 : std::vector<Int> s(fs.begin(), fs.end());
1346 676 : dgridsd(xyPositions.getStorage(del),
1347 : datStorage,
1348 338 : &s[0],
1349 338 : &s[1],
1350 338 : flags.getStorage(del),
1351 338 : rowFlags.getStorage(del),
1352 338 : &s[2],
1353 : &row,
1354 338 : griddedData.getStorage(del),
1355 : &nx,
1356 : &ny,
1357 : &npol,
1358 : &nchan,
1359 : &convSupport,
1360 : &convSampling,
1361 : convFunc.getStorage(del),
1362 : chanMap.getStorage(del),
1363 : polMap.getStorage(del));
1364 :
1365 338 : data.putStorage(datStorage, isCopy);
1366 : }
1367 338 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1368 : }
1369 :
1370 : // Helper functions for SDGrid::makeImage
1371 : namespace {
1372 : inline
1373 4322 : void setupVisBufferForFTMachineType(FTMachine::Type type, VisBuffer& vb) {
1374 4322 : switch(type) {
1375 0 : case FTMachine::RESIDUAL:
1376 0 : vb.visCube() = vb.correctedVisCube();
1377 0 : vb.visCube() -= vb.modelVisCube();
1378 0 : break;
1379 0 : case FTMachine::PSF:
1380 0 : vb.visCube() = Complex(1.0,0.0);
1381 0 : break;
1382 2161 : case FTMachine::COVERAGE:
1383 2161 : vb.visCube() = Complex(1.0);
1384 2161 : break;
1385 2161 : default:
1386 2161 : break;
1387 : }
1388 4322 : }
1389 :
1390 : inline
1391 258 : void getParamsForFTMachineType(const ROVisibilityIterator& vi, FTMachine::Type in_type,
1392 : casacore::Bool& out_dopsf, FTMachine::Type& out_type) {
1393 :
1394 : // Tune input type of FTMachine
1395 258 : auto haveCorrectedData = not (vi.msColumns().correctedData().isNull());
1396 258 : auto tunedType =
1397 258 : ((in_type == FTMachine::CORRECTED) && (not haveCorrectedData)) ?
1398 : FTMachine::OBSERVED : in_type;
1399 :
1400 : // Compute output parameters
1401 258 : switch(tunedType) {
1402 65 : case FTMachine::RESIDUAL:
1403 : case FTMachine::MODEL:
1404 : case FTMachine::CORRECTED:
1405 65 : out_dopsf = false;
1406 65 : out_type = tunedType;
1407 65 : break;
1408 129 : case FTMachine::PSF:
1409 : case FTMachine::COVERAGE:
1410 129 : out_dopsf = true;
1411 129 : out_type = tunedType;
1412 129 : break;
1413 64 : default:
1414 64 : out_dopsf = false;
1415 64 : out_type = FTMachine::OBSERVED;
1416 64 : break;
1417 : }
1418 258 : }
1419 :
1420 1186 : void abortOnPolFrameChange(const StokesImageUtil::PolRep refPolRep, const String & refMsName, ROVisibilityIterator &vi) {
1421 1186 : auto vb = vi.getVisBuffer();
1422 1186 : const auto polRep = (vb->polFrame() == MSIter::Linear) ?
1423 1186 : StokesImageUtil::LINEAR : StokesImageUtil::CIRCULAR;
1424 1186 : const auto polFrameChanged = (polRep != refPolRep);
1425 1186 : if (polFrameChanged) {
1426 : // Time of polarization change
1427 0 : constexpr auto timeColEnum = MS::PredefinedColumns::TIME;
1428 0 : const auto & timeColName = MS::columnName(timeColEnum);
1429 0 : const auto timeVbFirstRow = vb->time()[0];
1430 :
1431 0 : ScalarMeasColumn<MEpoch> timeMeasCol(vi.ms(),timeColName);
1432 0 : const auto & timeMeasRef = timeMeasCol.getMeasRef();
1433 :
1434 0 : const auto & timeUnit = MS::columnUnit(timeColEnum);
1435 :
1436 0 : MEpoch polFrameChangeEpoch(Quantity(timeVbFirstRow,timeUnit),
1437 0 : timeMeasRef);
1438 0 : MVTime polFrameChangeTime(polFrameChangeEpoch.getValue());
1439 :
1440 : // Error message
1441 0 : MVTime::Format fmt(MVTime::YMD | MVTime::USE_SPACE,10);
1442 0 : constexpr auto nl = '\n';
1443 0 : stringstream msg;
1444 : msg << "Polarization Frame changed ! " << nl
1445 0 : << setw(9) << right << "from: " << setw(8) << left << StokesImageUtil::toString(refPolRep)
1446 0 : << " in: " << refMsName << nl
1447 0 : << setw(9) << right << "to: " << setw(8) << left << StokesImageUtil::toString(polRep)
1448 0 : << " at: " << MVTime::Format(fmt) << polFrameChangeTime
1449 0 : << " (" << fixed << setprecision(6) << polFrameChangeTime.second() << " s)"
1450 0 : << " in: " << vb->msName();
1451 :
1452 0 : LogOrigin logOrigin("","abortOnPolFrameChange()");
1453 0 : LogIO logger(logOrigin);
1454 0 : logger << msg.str() << LogIO::SEVERE << LogIO::POST;
1455 :
1456 : // Abort
1457 0 : logger.sourceLocation(WHERE);
1458 0 : logger << "Polarization Frame must not change" << LogIO::EXCEPTION;
1459 : };
1460 1186 : }
1461 :
1462 : } // SDGrid::makeImage helper functions namespace
1463 :
1464 1186 : void SDGrid::nextChunk(ROVisibilityIterator &vi) {
1465 : #if defined(SDGRID_PERFS)
1466 2372 : StartStop trigger(cNextChunk);
1467 : #endif
1468 1186 : vi.nextChunk();
1469 1186 : }
1470 :
1471 : // Make a plain straightforward honest-to-FSM image. This returns
1472 : // a complex image, without conversion to Stokes. The representation
1473 : // is that required for the visibilities.
1474 : //----------------------------------------------------------------------
1475 228 : void SDGrid::makeImage(FTMachine::Type inType,
1476 : ROVisibilityIterator& vi,
1477 : ImageInterface<Complex>& theImage,
1478 : Matrix<Float>& weight) {
1479 :
1480 : // Attach visibility buffer (VisBuffer) to visibility iterator (VisibilityIterator)
1481 456 : VisBuffer vb(vi);
1482 :
1483 : // Set the Polarization Representation
1484 : // of image's Stokes Coordinate Axis
1485 : // based on first data in first MS
1486 228 : vi.origin();
1487 228 : const auto imgPolRep = (vb.polFrame() == MSIter::Linear) ?
1488 228 : StokesImageUtil::LINEAR : StokesImageUtil::CIRCULAR;
1489 228 : StokesImageUtil::changeCStokesRep(theImage, imgPolRep);
1490 456 : const auto firstMsName = vb.msName();
1491 :
1492 228 : initializeToSky(theImage, weight, vb);
1493 :
1494 : // Setup SDGrid Cache Manager
1495 228 : const auto onDuty = cacheIsEnabled;
1496 228 : const auto accessMode = cache.isEmpty() ? Cache::AccessMode::WRITE
1497 228 : : Cache::AccessMode::READ;
1498 456 : CacheManager cacheManager(cache, onDuty, accessMode);
1499 :
1500 : // Loop over the visibilities, putting VisBuffers
1501 1414 : for (vi.originChunks(); vi.moreChunks(); nextChunk(vi)) {
1502 1186 : abortOnPolFrameChange(imgPolRep, firstMsName, vi);
1503 : FTMachine::Type actualType;
1504 : Bool doPSF;
1505 1186 : if (vi.newMS()) { // Note: the first MS is a new MS
1506 258 : getParamsForFTMachineType(vi, inType, doPSF, actualType);
1507 258 : handleNewMs(vi, theImage);
1508 : }
1509 5508 : for (vi.origin(); vi.more(); vi++) {
1510 4322 : setupVisBufferForFTMachineType(actualType, vb);
1511 4322 : constexpr Int allVbRows = -1;
1512 4322 : put(vb, allVbRows, doPSF, actualType);
1513 : }
1514 : }
1515 228 : finalizeToSky();
1516 :
1517 : // Normalize by dividing out weights, etc.
1518 228 : auto doNormalize = (inType != FTMachine::COVERAGE);
1519 228 : getImage(weight, doNormalize);
1520 :
1521 : // Warning message
1522 228 : if (allEQ(weight, 0.0f)) {
1523 6 : LogIO logger(LogOrigin(name(),"makeImage"));
1524 : logger << LogIO::SEVERE
1525 : << "No useful data in SDGrid: all weights are zero"
1526 2 : << LogIO::POST;
1527 : }
1528 228 : }
1529 :
1530 : // Finalize : optionally normalize by weight image
1531 228 : ImageInterface<Complex>& SDGrid::getImage(Matrix<Float>& weights,
1532 : Bool normalize)
1533 : {
1534 228 : AlwaysAssert(lattice, AipsError);
1535 228 : AlwaysAssert(image, AipsError);
1536 :
1537 228 : logIO() << LogOrigin("SDGrid", "getImage") << LogIO::NORMAL;
1538 :
1539 : // execute minmax clipping
1540 228 : clipMinMax();
1541 :
1542 228 : weights.resize(sumWeight.shape());
1543 :
1544 228 : convertArray(weights,sumWeight);
1545 :
1546 : // If the weights are all zero then we cannot normalize
1547 : // otherwise we don't care.
1548 : ///////////////////////
1549 : /*{
1550 : PagedImage<Float> thisScreen(lattice->shape(), image->coordinates(), "TheData");
1551 : LatticeExpr<Float> le(abs(*lattice));
1552 : thisScreen.copyData(le);
1553 : PagedImage<Float> thisScreen2(lattice->shape(), image->coordinates(), "TheWeight");
1554 : LatticeExpr<Float> le2(abs(*wLattice));
1555 : thisScreen2.copyData(le2);
1556 : }*/
1557 : /////////////////////
1558 :
1559 228 : if(normalize) {
1560 114 : if(max(weights)==0.0) {
1561 : //logIO() << LogIO::WARN << "No useful data in SDGrid: weights all zero"
1562 : // << LogIO::POST;
1563 : //image->set(Complex(0.0));
1564 1 : return *image;
1565 : }
1566 : //Timer tim;
1567 : //tim.mark();
1568 : //lattice->copyData((LatticeExpr<Complex>)(iif((*wLattice<=minWeight_p), 0.0,
1569 : // (*lattice)/(*wLattice))));
1570 : //As we are not using disk based lattices anymore the above uses too much memory for
1571 : // ArrayLattice...it does not do a real inplace math but rather mutiple copies
1572 : // seem to be involved thus the less elegant but faster and less memory hog loop
1573 : // below.
1574 :
1575 226 : IPosition pos(4);
1576 7326 : for (Int chan=0; chan < nchan; ++chan){
1577 7213 : pos[3]=chan;
1578 14467 : for( Int pol=0; pol < npol; ++ pol){
1579 7254 : pos[2]=pol;
1580 455875 : for (Int y=0; y < ny ; ++y){
1581 448621 : pos[1]=y;
1582 34433263 : for (Int x=0; x < nx; ++x){
1583 33984642 : pos[0]=x;
1584 33984642 : Float wgt=wGriddedData(pos);
1585 33984642 : if(wgt > minWeight_p)
1586 26113193 : griddedData(pos)=griddedData(pos)/wgt;
1587 : else
1588 7871449 : griddedData(pos)=0.0;
1589 : }
1590 : }
1591 : }
1592 : }
1593 : //tim.show("After normalizing");
1594 : }
1595 :
1596 : //if(!isTiled)
1597 : {
1598 : // Now find the SubLattice corresponding to the image
1599 454 : IPosition gridShape(4, nx, ny, npol, nchan);
1600 454 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
1601 681 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1602 454 : IPosition stride(4, 1);
1603 681 : IPosition trc(blc+image->shape()-stride);
1604 454 : LCBox gridBox(blc, trc, gridShape);
1605 681 : SubLattice<Complex> gridSub(*arrayLattice, gridBox);
1606 :
1607 : // Do the copy
1608 227 : image->copyData(gridSub);
1609 : }
1610 227 : return *image;
1611 : }
1612 :
1613 : // Return weights image
1614 0 : void SDGrid::getWeightImage(ImageInterface<Float>& weightImage, Matrix<Float>& weights)
1615 : {
1616 0 : AlwaysAssert(lattice, AipsError);
1617 0 : AlwaysAssert(image, AipsError);
1618 :
1619 0 : logIO() << LogOrigin("SDGrid", "getWeightImage") << LogIO::NORMAL;
1620 :
1621 0 : weights.resize(sumWeight.shape());
1622 0 : convertArray(weights,sumWeight);
1623 :
1624 0 : weightImage.copyData(*wArrayLattice);
1625 0 : }
1626 :
1627 476 : void SDGrid::ok() {
1628 476 : AlwaysAssert(image, AipsError);
1629 476 : }
1630 :
1631 : // Get the index into the pointing table for this time. Note that the
1632 : // in the pointing table, TIME specifies the beginning of the spanned
1633 : // time range, whereas for the main table, TIME is the centroid.
1634 : // Note that the behavior for multiple matches is not specified! i.e.
1635 : // if there are multiple matches, the index returned depends on the
1636 : // history of previous matches. It is deterministic but not obvious.
1637 : // One could cure this by searching but it would be considerably
1638 : // costlier.
1639 179093 : Int SDGrid::getIndex(const MSPointingColumns& mspc, const Double& time,
1640 : const Double& interval, const Int& antid) {
1641 : //Int start=lastIndex_p;
1642 179093 : Int start=lastIndexPerAnt_p[antid];
1643 179093 : Double tol=interval < 0.0 ? time*C::dbl_epsilon : interval/2.0;
1644 : // Search forwards
1645 179093 : Int nrows=mspc.time().nrow();
1646 402416 : for (Int i=start;i<nrows;i++) {
1647 : // Check for ANTENNA_ID. When antid<0 (default) ANTENNA_ID in POINTING table < 0 is ignored.
1648 : // MS v2 definition indicates ANTENNA_ID in POINING >= 0.
1649 402398 : if (antid >= 0 && mspc.antennaId()(i) != antid) {
1650 0 : continue;
1651 : }
1652 :
1653 402398 : Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
1654 : // If the interval in the pointing table is negative, use the last
1655 : // entry. Note that this may be invalid (-1) but in that case
1656 : // the calling routine will generate an error
1657 402398 : Double mspc_interval = mspc.interval()(i);
1658 :
1659 402398 : if(mspc_interval<0.0) {
1660 : //return lastIndex_p;
1661 0 : return lastIndexPerAnt_p[antid];
1662 : }
1663 : // Pointing table interval is specified so we have to do a match
1664 : else {
1665 : // Is the midpoint of this pointing table entry within the specified
1666 : // tolerance of the main table entry?
1667 402398 : if(abs(midpoint-time) <= (mspc_interval/2.0+tol)) {
1668 : //lastIndex_p=i;
1669 179075 : lastIndexPerAnt_p[antid]=i;
1670 179075 : return i;
1671 : }
1672 : }
1673 : }
1674 : // Search backwards
1675 69174 : for (Int i=start;i>=0;i--) {
1676 69174 : if (antid >= 0 && mspc.antennaId()(i) != antid) {
1677 0 : continue;
1678 : }
1679 69174 : Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
1680 69174 : Double mspc_interval = mspc.interval()(i);
1681 69174 : if(mspc_interval<0.0) {
1682 : //return lastIndex_p;
1683 0 : return lastIndexPerAnt_p[antid];
1684 : }
1685 : // Pointing table interval is specified so we have to do a match
1686 : else {
1687 : // Is the midpoint of this pointing table entry within the specified
1688 : // tolerance of the main table entry?
1689 69174 : if(abs(midpoint-time) <= (mspc_interval/2.0+tol)) {
1690 : //lastIndex_p=i;
1691 18 : lastIndexPerAnt_p[antid]=i;
1692 18 : return i;
1693 : }
1694 : }
1695 : }
1696 : // No match!
1697 0 : return -1;
1698 : }
1699 :
1700 355086 : Bool SDGrid::getXYPos(const VisBuffer& vb, Int row) {
1701 :
1702 : // Cache control
1703 355086 : if (cacheIsEnabled and cache.isReadable()) {
1704 176039 : cache.loadRowPixel();
1705 176039 : return rowPixel.isValid;
1706 : }
1707 179047 : const auto onDuty = cacheIsEnabled and cache.isWriteable();
1708 358094 : CacheWriter cacheWriter(cache, onDuty);
1709 :
1710 : // Until we manage to compute a valid one ...
1711 179047 : rowPixel.isValid = false;
1712 :
1713 : // Select the POINTING table (columns) we'll work with.
1714 179047 : const auto haveConvertedColumn = ramPointingTable.nrow() > 0;
1715 179047 : const MSPointingColumns& pointingColumns = haveConvertedColumn ?
1716 0 : *ramPointingColumnsPtr
1717 179047 : : vb.msColumns().pointing();
1718 :
1719 179047 : const auto nPointings = pointingColumns.nrow();
1720 179047 : const auto havePointings = (nPointings >= 1);
1721 :
1722 : // We'll need to call these many times, so let's call them once for good
1723 179047 : const auto rowTime = vb.time()(row);
1724 179047 : const auto rowTimeInterval = vb.timeInterval()(row);
1725 179047 : const auto rowAntenna1 = vb.antenna1()(row);
1726 :
1727 : // 1. Try to find the index of a pointing recorded:
1728 : // - for the antenna of specified row,
1729 : // - at a time close enough to the time at which
1730 : // data of specified row was taken using that antenna
1731 179047 : constexpr Int invalidIndex = -1;
1732 179047 : Int pointingIndex = invalidIndex;
1733 179047 : if (havePointings) {
1734 : #if defined(SDGRID_PERFS)
1735 179047 : StartStop trigger(cSearchValidPointing);
1736 : #endif
1737 : // if (vb.newMS()) vb.newMS does not work well using msid
1738 : // Note about above comment:
1739 : // - vb.newMS probably works well
1740 : // - but if the calling code is iterating over the rows of a subchunk
1741 : // vb.newMS returns true for all rows belonging to the first subchunk
1742 : // of the first chunk of a new MS.
1743 :
1744 : // ???
1745 : // What if vb changed since we were last called ?
1746 : // What if the calling code calls put and get, with different VisBuffers ?
1747 179047 : if (vb.msId() != msId_p) {
1748 139 : lastIndex_p = 0; // no longer used
1749 139 : if (lastIndexPerAnt_p.nelements() < (size_t)vb.numberAnt()) {
1750 96 : lastIndexPerAnt_p.resize(vb.numberAnt());
1751 : }
1752 139 : lastIndexPerAnt_p = 0;
1753 139 : msId_p = vb.msId();
1754 139 : lastAntID_p = -1;
1755 : }
1756 :
1757 : // Try to locate a pointing verifying for specified row:
1758 : // | POINTING.TIME - MAIN.TIME | <= 0.5*(POINTING.INTERVAL + tolerance)
1759 :
1760 : // Try first using a tiny tolerance
1761 179047 : constexpr Double useTinyTolerance = -1.0;
1762 179047 : pointingIndex = getIndex(pointingColumns, rowTime, useTinyTolerance , rowAntenna1);
1763 :
1764 179047 : auto foundPointing = (pointingIndex >= 0);
1765 179047 : if (not foundPointing) {
1766 : // Try again using tolerance = MAIN.INTERVAL
1767 0 : pointingIndex = getIndex(pointingColumns, rowTime, rowTimeInterval, rowAntenna1);
1768 0 : foundPointing = (pointingIndex >= 0);
1769 : }
1770 :
1771 : // Making the implicit type conversion explicit.
1772 : // Conversion is safe because it occurs only when pointingIndex >= 0.
1773 179047 : const auto foundValidPointing = (foundPointing and (static_cast<rownr_t>(pointingIndex) < nPointings));
1774 :
1775 179047 : if (not foundValidPointing) {
1776 0 : ostringstream o;
1777 0 : o << "Failed to find pointing information for time "
1778 0 : << MVTime(rowTime/86400.0) << " : Omitting this point";
1779 :
1780 0 : logIO_p << LogIO::DEBUGGING << String(o) << LogIO::POST;
1781 :
1782 0 : return rowPixel.isValid;
1783 : }
1784 : }
1785 :
1786 : // 2. At this stage we have:
1787 : // * either no pointings and an invalid pointingIndex
1788 : // * or pointings and a valid pointingIndex.
1789 : // Decide now if we need to interpolate antenna's pointing direction
1790 : // at data-taking time:
1791 : // we'll do so when data is sampled faster than pointings are recorded
1792 179047 : Bool needInterpolation = False;
1793 179047 : if (havePointings) {
1794 179047 : const auto pointingInterval = pointingColumns.interval()(pointingIndex);
1795 179047 : if (rowTimeInterval < pointingInterval) needInterpolation = True;
1796 : }
1797 179047 : const auto mustInterpolate = havePointings && needInterpolation;
1798 :
1799 : // 3. Create interpolator if needed
1800 179047 : if (mustInterpolate) {
1801 5000 : if (not isSplineInterpolationReady) {
1802 : #if defined(SDGRID_PERFS)
1803 3 : StartStop trigger(cComputeSplines);
1804 : #endif
1805 : const auto nAntennas = static_cast<size_t>(
1806 3 : vb.msColumns().antenna().nrow()
1807 : );
1808 3 : interpolator = new SDPosInterpolator(
1809 : pointingColumns,
1810 3 : pointingDirCol_p,
1811 : nAntennas
1812 3 : );
1813 3 : isSplineInterpolationReady = true;
1814 : } else {
1815 : // We have an interpolator. Re-use it if possible.
1816 4997 : const auto canReuseInterpolator = interpolator->inTimeRange(rowTime, rowAntenna1);
1817 4997 : if (not canReuseInterpolator) {
1818 : // setup spline interpolator for the current dataset (CAS-11261, 2018/5/22 WK)
1819 : // delete and re-create it
1820 2 : delete interpolator;
1821 2 : interpolator = 0;
1822 : const auto nAntennas = static_cast<size_t>(
1823 2 : vb.msColumns().antenna().nrow()
1824 : );
1825 2 : interpolator = new SDPosInterpolator(
1826 : pointingColumns,
1827 2 : pointingDirCol_p,
1828 : nAntennas
1829 2 : );
1830 : }
1831 : }
1832 : }
1833 :
1834 : // 4. Create the direction conversion machine if needed
1835 :
1836 358094 : if ( pointingDirCol_p == "SOURCE_OFFSET" or
1837 179047 : pointingDirCol_p == "POINTING_OFFSET" ) {
1838 : // It makes no sense to track in offset coordinates...
1839 : // hopefully the user sets the image coords right
1840 0 : fixMovingSource_p = false;
1841 : }
1842 :
1843 179047 : const auto needDirectionConverter = (
1844 179047 : not havePointings or not haveConvertedColumn or fixMovingSource_p
1845 : );
1846 :
1847 179047 : if (needDirectionConverter) {
1848 179047 : if (not pointingToImage) {
1849 : // Set the frame
1850 : const auto & rowAntenna1Position =
1851 248 : vb.msColumns().antenna().positionMeas()(rowAntenna1);
1852 : // set dummy time stamp 1 day before rowTime
1853 372 : const MEpoch dummyEpoch(Quantity(rowTime - 86400.0, "s"));
1854 124 : mFrame_p = MeasFrame(dummyEpoch, rowAntenna1Position);
1855 :
1856 : // Remember antenna id for next call,
1857 : // which may be done using a different VisBuffer ...
1858 124 : lastAntID_p = rowAntenna1;
1859 :
1860 : // Compute the "model" required to setup the conversion machine
1861 124 : if (havePointings) {
1862 248 : worldPosMeas = mustInterpolate ?
1863 : directionMeas(pointingColumns, pointingIndex, rowTime)
1864 124 : : directionMeas(pointingColumns, pointingIndex);
1865 : } else {
1866 : // Without pointings, this sets the direction to the phase center
1867 0 : worldPosMeas = vb.direction1()(row);
1868 : }
1869 :
1870 : // Make a machine to convert from the worldPosMeas to the output
1871 : // Direction Measure type for the relevant frame
1872 248 : MDirection::Ref outRef(directionCoord.directionType(), mFrame_p);
1873 124 : pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
1874 124 : if (not pointingToImage) {
1875 0 : logIO_p << "Cannot make direction conversion machine" << LogIO::EXCEPTION;
1876 : }
1877 : // Perform 1 dummy direction conversion to clear values
1878 : // cached in static variables of casacore functions like MeasTable::dUT1
1879 124 : MDirection _dir_tmp = (*pointingToImage)();
1880 : }
1881 : }
1882 :
1883 : // 5. Update the frame holding the measurements for this row
1884 537141 : const MEpoch rowEpoch(Quantity(rowTime, "s"));
1885 : // ---- Always reset the epoch
1886 : {
1887 : #if defined(SDGRID_PERFS)
1888 358094 : StartStop trigger(cResetFrame);
1889 : #endif
1890 179047 : mFrame_p.resetEpoch(rowEpoch);
1891 : }
1892 : // ---- Reset antenna position only if antenna changed since we were last called
1893 179047 : if (lastAntID_p != rowAntenna1) {
1894 : // Debug messages
1895 15 : if (lastAntID_p == -1) {
1896 : // antenna ID is unset
1897 15 : logIO_p << LogIO::DEBUGGING
1898 : << "updating antenna position for conversion: new MS ID " << msId_p
1899 15 : << ", antenna ID " << rowAntenna1 << LogIO::POST;
1900 : } else {
1901 0 : logIO_p << LogIO::DEBUGGING
1902 : << "updating antenna position for conversion: MS ID " << msId_p
1903 : << ", last antenna ID " << lastAntID_p
1904 0 : << ", new antenna ID " << rowAntenna1 << LogIO::POST;
1905 : }
1906 : const auto & rowAntenna1Position =
1907 15 : vb.msColumns().antenna().positionMeas()(rowAntenna1);
1908 : {
1909 : #if defined(SDGRID_PERFS)
1910 30 : StartStop trigger(cResetFrame);
1911 : #endif
1912 15 : mFrame_p.resetPosition(rowAntenna1Position);
1913 : }
1914 : // Remember antenna id for next call,
1915 : // which may be done using a different VisBuffer ...
1916 15 : lastAntID_p = rowAntenna1;
1917 : }
1918 :
1919 : // 6. Compute user-specified column direction at data-taking time,
1920 : // in image's direction reference frame
1921 179047 : if (havePointings) {
1922 179047 : if (mustInterpolate) {
1923 : #if defined(SDGRID_PERFS)
1924 5000 : cInterpolateDirection.start();
1925 : #endif
1926 : const auto interpolatedDirection =
1927 10000 : directionMeas(pointingColumns, pointingIndex, rowTime);
1928 : #if defined(SDGRID_PERFS)
1929 5000 : cInterpolateDirection.stop();
1930 5000 : cConvertDirection.start();
1931 : #endif
1932 : worldPosMeas = haveConvertedColumn ?
1933 : interpolatedDirection
1934 5000 : : (*pointingToImage)(interpolatedDirection);
1935 : #if defined(SDGRID_PERFS)
1936 5000 : cConvertDirection.stop();
1937 : #endif
1938 :
1939 : // Debug stuff
1940 : //Vector<Double> newdirv = newdir.getAngle("rad").getValue();
1941 : //cerr<<"dir0="<<newdirv(0)<<endl;
1942 :
1943 : //fprintf(pfile,"%.8f %.8f \n", newdirv(0), newdirv(1));
1944 : //printf("%lf %lf \n", newdirv(0), newdirv(1));
1945 : } else {
1946 348094 : const auto columnDirection = directionMeas(pointingColumns, pointingIndex);
1947 : worldPosMeas = haveConvertedColumn ?
1948 : columnDirection
1949 174047 : : (*pointingToImage)(columnDirection);
1950 : }
1951 : } else {
1952 : // Without pointings, this converts the direction of the phase center
1953 0 : worldPosMeas = (*pointingToImage)(vb.direction1()(row));
1954 : }
1955 :
1956 : // 7. Convert world position coordinates to image pixel coordinates
1957 : {
1958 : #if defined(SDGRID_PERFS)
1959 179047 : StartStop trigger(cComputeDirectionPixel);
1960 : #endif
1961 :
1962 179047 : rowPixel.isValid = directionCoord.toPixel(xyPos, worldPosMeas);
1963 : }
1964 :
1965 179047 : if (not rowPixel.isValid) {
1966 0 : logIO_p << "Failed to find a pixel for pointing direction of "
1967 0 : << MVTime(worldPosMeas.getValue().getLong("rad")).string(MVTime::TIME)
1968 0 : << ", " << MVAngle(worldPosMeas.getValue().getLat("rad")).string(MVAngle::ANGLE)
1969 0 : << LogIO::WARN << LogIO::POST;
1970 0 : return rowPixel.isValid;
1971 : }
1972 :
1973 : // 8. Handle moving sources
1974 179047 : if (fixMovingSource_p) {
1975 : #if defined(SDGRID_PERFS)
1976 19272 : StartStop trigger(cHandleMovingSource);
1977 : #endif
1978 9636 : if (xyPosMovingOrig_p.nelements() < 2) {
1979 1 : directionCoord.toPixel(xyPosMovingOrig_p, firstMovingDir_p);
1980 : }
1981 : //via HADEC or AZEL for parallax of near sources
1982 19272 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
1983 19272 : MDirection tmphadec = MDirection::Convert(movingDir_p, outref1)();
1984 19272 : MDirection actSourceDir = (*pointingToImage)(tmphadec);
1985 9636 : Vector<Double> actPix;
1986 9636 : directionCoord.toPixel(actPix, actSourceDir);
1987 :
1988 : //cout << row << " scan " << vb.scan()(row) << "xyPos " << xyPos
1989 : // << " xyposmovorig " << xyPosMovingOrig_p << " actPix " << actPix << endl;
1990 :
1991 9636 : xyPos = xyPos + xyPosMovingOrig_p - actPix;
1992 : }
1993 :
1994 179047 : return rowPixel.isValid;
1995 : }
1996 :
1997 174260 : MDirection SDGrid::directionMeas(const MSPointingColumns& mspc, const Int& index){
1998 174260 : if (pointingDirCol_p == "TARGET") {
1999 0 : return mspc.targetMeas(index);
2000 174260 : } else if (pointingDirCol_p == "POINTING_OFFSET") {
2001 0 : if (!mspc.pointingOffsetMeasCol().isNull()) {
2002 0 : return mspc.pointingOffsetMeas(index);
2003 : }
2004 0 : cerr << "No PONTING_OFFSET column in POINTING table" << endl;
2005 174260 : } else if (pointingDirCol_p == "SOURCE_OFFSET") {
2006 0 : if (!mspc.sourceOffsetMeasCol().isNull()) {
2007 0 : return mspc.sourceOffsetMeas(index);
2008 : }
2009 0 : cerr << "No SOURCE_OFFSET column in POINTING table" << endl;
2010 174260 : } else if (pointingDirCol_p == "ENCODER") {
2011 0 : if (!mspc.encoderMeas().isNull()) {
2012 0 : return mspc.encoderMeas()(index);
2013 : }
2014 0 : cerr << "No ENCODER column in POINTING table" << endl;
2015 : }
2016 :
2017 : //default return this
2018 174260 : return mspc.directionMeas(index);
2019 : }
2020 :
2021 : // for the cases, interpolation of the pointing direction requires
2022 : // when data sampling rate higher than the pointing data recording
2023 : // (e.g. fast OTF)
2024 5003 : MDirection SDGrid::directionMeas(const MSPointingColumns& mspc, const Int& index,
2025 : const Double& time){
2026 : //spline interpolation
2027 5003 : if (isSplineInterpolationReady) {
2028 5003 : Int antid = mspc.antennaId()(index);
2029 5003 : if (interpolator->doSplineInterpolation(antid)) {
2030 5003 : return interpolator->interpolateDirectionMeasSpline(mspc, time, index, antid);
2031 : }
2032 : }
2033 :
2034 : //linear interpolation (as done before CAS-7911)
2035 : Int index1, index2;
2036 0 : if (time < mspc.time()(index)) {
2037 0 : if (index > 0) {
2038 0 : index1 = index-1;
2039 0 : index2 = index;
2040 0 : } else if (index == 0) {
2041 0 : index1 = index;
2042 0 : index2 = index+1;
2043 : }
2044 : } else {
2045 0 : if (index < Int(mspc.nrow()-1)) {
2046 0 : index1 = index;
2047 0 : index2 = index+1;
2048 0 : } else if (index == Int(mspc.nrow()-1) || (mspc.time()(index)-mspc.time()(index+1))>2*mspc.interval()(index)) {
2049 0 : index1 = index-1;
2050 0 : index2 = index;
2051 : }
2052 : }
2053 0 : return interpolateDirectionMeas(mspc, time, index, index1, index2);
2054 : }
2055 :
2056 0 : MDirection SDGrid::interpolateDirectionMeas(const MSPointingColumns& mspc,
2057 : const Double& time,
2058 : const Int& index,
2059 : const Int& indx1,
2060 : const Int& indx2){
2061 0 : Vector<Double> dir1,dir2;
2062 0 : Vector<Double> newdir(2),scanRate(2);
2063 : Double dLon, dLat;
2064 : Double ftime,ftime2,ftime1,dtime;
2065 0 : MDirection newDirMeas;
2066 0 : MDirection::Ref rf;
2067 : Bool isfirstRefPt;
2068 :
2069 0 : if (indx1 == index) {
2070 0 : isfirstRefPt = true;
2071 : } else {
2072 0 : isfirstRefPt = false;
2073 : }
2074 :
2075 0 : if (pointingDirCol_p == "TARGET") {
2076 0 : dir1 = mspc.targetMeas(indx1).getAngle("rad").getValue();
2077 0 : dir2 = mspc.targetMeas(indx2).getAngle("rad").getValue();
2078 0 : } else if (pointingDirCol_p == "POINTING_OFFSET") {
2079 0 : if (!mspc.pointingOffsetMeasCol().isNull()) {
2080 0 : dir1 = mspc.pointingOffsetMeas(indx1).getAngle("rad").getValue();
2081 0 : dir2 = mspc.pointingOffsetMeas(indx2).getAngle("rad").getValue();
2082 : } else {
2083 0 : cerr << "No PONTING_OFFSET column in POINTING table" << endl;
2084 : }
2085 0 : } else if (pointingDirCol_p == "SOURCE_OFFSET") {
2086 0 : if (!mspc.sourceOffsetMeasCol().isNull()) {
2087 0 : dir1 = mspc.sourceOffsetMeas(indx1).getAngle("rad").getValue();
2088 0 : dir2 = mspc.sourceOffsetMeas(indx2).getAngle("rad").getValue();
2089 : } else {
2090 0 : cerr << "No SOURCE_OFFSET column in POINTING table" << endl;
2091 : }
2092 0 : } else if (pointingDirCol_p == "ENCODER") {
2093 0 : if (!mspc.encoderMeas().isNull()) {
2094 0 : dir1 = mspc.encoderMeas()(indx1).getAngle("rad").getValue();
2095 0 : dir2 = mspc.encoderMeas()(indx2).getAngle("rad").getValue();
2096 : } else {
2097 0 : cerr << "No ENCODER column in POINTING table" << endl;
2098 : }
2099 : } else {
2100 0 : dir1 = mspc.directionMeas(indx1).getAngle("rad").getValue();
2101 0 : dir2 = mspc.directionMeas(indx2).getAngle("rad").getValue();
2102 : }
2103 :
2104 0 : dLon = dir2(0) - dir1(0);
2105 0 : dLat = dir2(1) - dir1(1);
2106 0 : ftime = floor(mspc.time()(indx1));
2107 0 : ftime2 = mspc.time()(indx2) - ftime;
2108 0 : ftime1 = mspc.time()(indx1) - ftime;
2109 0 : dtime = ftime2 - ftime1;
2110 0 : scanRate(0) = dLon/dtime;
2111 0 : scanRate(1) = dLat/dtime;
2112 : //scanRate(0) = dir2(0)/dtime-dir1(0)/dtime;
2113 : //scanRate(1) = dir2(1)/dtime-dir1(1)/dtime;
2114 : //Double delT = mspc.time()(index)-time;
2115 : //cerr<<"index="<<index<<" dLat="<<dLat<<" dtime="<<dtime<<" delT="<< delT<<endl;
2116 : //cerr<<"deldirlat="<<scanRate(1)*fabs(delT)<<endl;
2117 0 : if (isfirstRefPt) {
2118 0 : newdir(0) = dir1(0)+scanRate(0)*fabs(mspc.time()(index)-time);
2119 0 : newdir(1) = dir1(1)+scanRate(1)*fabs(mspc.time()(index)-time);
2120 0 : rf = mspc.directionMeas(indx1).getRef();
2121 : } else {
2122 0 : newdir(0) = dir2(0)-scanRate(0)*fabs(mspc.time()(index)-time);
2123 0 : newdir(1) = dir2(1)-scanRate(1)*fabs(mspc.time()(index)-time);
2124 0 : rf = mspc.directionMeas(indx2).getRef();
2125 : }
2126 : //default return this
2127 0 : Quantity rDirLon(newdir(0), "rad");
2128 0 : Quantity rDirLat(newdir(1), "rad");
2129 0 : newDirMeas = MDirection(rDirLon, rDirLat, rf);
2130 : //cerr<<"newDirMeas rf="<<newDirMeas.getRefString()<<endl;
2131 : //return mspc.directionMeas(index);
2132 0 : return newDirMeas;
2133 : }
2134 :
2135 4322 : void SDGrid::pickWeights(const VisBuffer& vb, Matrix<Float>& weight){
2136 : #if defined(SDGRID_PERFS)
2137 8644 : StartStop trigger(cPickWeights);
2138 : #endif
2139 : //break reference
2140 4322 : weight.resize();
2141 :
2142 4322 : if (useImagingWeight_p) {
2143 0 : weight.reference(vb.imagingWeight());
2144 : } else {
2145 8644 : const Cube<Float> weightspec(vb.weightSpectrum());
2146 4322 : weight.resize(vb.nChannel(), vb.nRow());
2147 :
2148 4322 : if (weightspec.nelements() == 0) {
2149 356400 : for (Int k = 0; k < vb.nRow(); ++k) {
2150 : //cerr << "nrow " << vb.nRow() << " " << weight.shape() << " " << weight.column(k).shape() << endl;
2151 352078 : weight.column(k).set(vb.weight()(k));
2152 : }
2153 : } else {
2154 0 : Int npol = weightspec.shape()(0);
2155 0 : if (npol == 1) {
2156 0 : for (Int k = 0; k < vb.nRow(); ++k) {
2157 0 : for (int chan = 0; chan < vb.nChannel(); ++chan) {
2158 0 : weight(chan, k)=weightspec(0, chan, k);
2159 : }
2160 : }
2161 : } else {
2162 0 : for (Int k = 0; k < vb.nRow(); ++k) {
2163 0 : for (int chan = 0; chan < vb.nChannel(); ++chan) {
2164 0 : weight(chan, k) = (weightspec(0, chan, k) + weightspec((npol-1), chan, k))/2.0f;
2165 : }
2166 : }
2167 : }
2168 : }
2169 : }
2170 4322 : }
2171 :
2172 228 : void SDGrid::clipMinMax() {
2173 228 : if (clipminmax_) {
2174 : Bool gmin_b, gmax_b, wmin_b, wmax_b, np_b;
2175 16 : const auto *gmin_p = gmin_.getStorage(gmin_b);
2176 16 : const auto *gmax_p = gmax_.getStorage(gmax_b);
2177 16 : const auto *wmin_p = wmin_.getStorage(wmin_b);
2178 16 : const auto *wmax_p = wmax_.getStorage(wmax_b);
2179 16 : const auto *np_p = npoints_.getStorage(np_b);
2180 :
2181 : Bool data_b, weight_b, sumw_b;
2182 16 : auto data_p = griddedData.getStorage(data_b);
2183 16 : auto weight_p = wGriddedData.getStorage(weight_b);
2184 16 : auto sumw_p = sumWeight.getStorage(sumw_b);
2185 :
2186 32 : auto arrayShape = griddedData.shape();
2187 16 : size_t num_xy = arrayShape.getFirst(2).product();
2188 16 : size_t num_polchan = arrayShape.getLast(2).product();
2189 160 : for (size_t i = 0; i < num_xy; ++i) {
2190 306 : for (size_t j = 0; j < num_polchan; ++j) {
2191 162 : auto k = i * num_polchan + j;
2192 162 : if (np_p[k] == 1) {
2193 2 : auto wt = wmin_p[k];
2194 2 : data_p[k] = wt * gmin_p[k];
2195 2 : weight_p[k] = wt;
2196 2 : sumw_p[j] += wt;
2197 160 : } else if (np_p[k] == 2) {
2198 1 : auto wt = wmin_p[k];
2199 1 : data_p[k] = wt * gmin_p[k];
2200 1 : weight_p[k] = wt;
2201 1 : sumw_p[j] += wt;
2202 1 : wt = wmax_p[k];
2203 1 : data_p[k] += wt * gmax_p[k];
2204 1 : weight_p[k] += wt;
2205 1 : sumw_p[j] += wt;
2206 : }
2207 : }
2208 : }
2209 :
2210 16 : wGriddedData.putStorage(weight_p, weight_b);
2211 16 : griddedData.putStorage(data_p, data_b);
2212 16 : sumWeight.putStorage(sumw_p, sumw_b);
2213 :
2214 16 : npoints_.freeStorage(np_p, np_b);
2215 16 : wmax_.freeStorage(wmax_p, wmax_b);
2216 16 : wmin_.freeStorage(wmin_p, wmin_b);
2217 16 : gmax_.freeStorage(gmax_p, gmax_b);
2218 16 : gmin_.freeStorage(gmin_p, gmin_b);
2219 : }
2220 228 : }
2221 :
2222 124 : SDGrid::MaskedPixelRef::MaskedPixelRef(Vector<Double>& xyIn, Bool isValidIn)
2223 : : xy {xyIn},
2224 124 : isValid {isValidIn}
2225 124 : {}
2226 :
2227 : SDGrid::MaskedPixelRef&
2228 0 : SDGrid::MaskedPixelRef::operator=(const SDGrid::MaskedPixelRef &other) {
2229 0 : xy = other.xy;
2230 0 : isValid = other.isValid;
2231 0 : return *this;
2232 : }
2233 :
2234 176039 : SDGrid::MaskedPixel::MaskedPixel(Double xIn, Double yIn, Bool isValidIn)
2235 : : x {xIn},
2236 : y {yIn},
2237 176039 : isValid {isValidIn}
2238 176039 : {}
2239 :
2240 124 : SDGrid::Cache::Cache(SDGrid &parent)
2241 : : sdgrid {parent},
2242 : isOpened {false},
2243 : accessMode {AccessMode::READ},
2244 : canRead {false},
2245 : canWrite {false},
2246 124 : inputPixel {sdgrid.rowPixel},
2247 124 : outputPixel {sdgrid.rowPixel}
2248 124 : {}
2249 :
2250 258 : const casacore::String& SDGrid::Cache::className() {
2251 258 : static casacore::String className_ = "SDGrid::Cache";
2252 258 : return className_;
2253 : }
2254 :
2255 0 : SDGrid::Cache& SDGrid::Cache::operator=(const Cache &other) {
2256 0 : sdgrid = other.sdgrid;
2257 0 : msCaches = other.msCaches;
2258 0 : isOpened = false;
2259 0 : canRead = false;
2260 0 : canWrite = false;
2261 : // inputPixel = sdgrid.rowPixel;
2262 0 : inputPixel.xy = sdgrid.rowPixel.xy;
2263 : //inputPixel.isValid = sdgrid.rowPixel.isValid;
2264 : // <=>
2265 0 : msPixels = nullptr;
2266 0 : outputPixel = other.outputPixel;
2267 0 : msCacheReadIterator = MsCaches::const_iterator();
2268 0 : pixelReadIterator = Pixels::const_iterator();
2269 0 : return *this;
2270 : }
2271 :
2272 : void
2273 228 : SDGrid::Cache::open(AccessMode accessModeIn) {
2274 228 : if (isOpened) {
2275 0 : LogIO logger(LogOrigin(className(), "open", WHERE));
2276 0 : logger << "BUG: Opened " << className() << " was re-opened before being closed." << LogIO::EXCEPTION;
2277 : }
2278 228 : isOpened = true;
2279 228 : accessMode = accessModeIn;
2280 228 : canRead = accessMode == AccessMode::READ;
2281 228 : canWrite = accessMode == AccessMode::WRITE;
2282 228 : if (outputPixel.xy.size() < 2) {
2283 114 : outputPixel.xy.resize(2);
2284 114 : outputPixel.isValid = false;
2285 : }
2286 228 : rewind();
2287 228 : }
2288 :
2289 : void
2290 228 : SDGrid::Cache::rewind() {
2291 : // Writer
2292 228 : msPixels = nullptr;
2293 :
2294 : // Reader
2295 228 : msCacheReadIterator = msCaches.cbegin();
2296 228 : pixelReadIterator = Pixels::const_iterator();
2297 228 : }
2298 :
2299 : void
2300 228 : SDGrid::Cache::close() {
2301 228 : if (not isOpened) {
2302 0 : LogIO logger(LogOrigin(className(), "close", WHERE));
2303 0 : logger << "BUG: Closed " << className() << " was re-closed before being opened." << LogIO::EXCEPTION;
2304 : }
2305 228 : if (isWriteable()) {
2306 : // Make sure we have 1 pixel per row
2307 243 : for (const auto & msCache : msCaches) {
2308 129 : if (not msCache.isConsistent()) {
2309 0 : LogIO logger(LogOrigin(className(), "close", WHERE));
2310 0 : const auto didOverflow = msCache.pixels.size() > msCache.nRows;
2311 0 : const auto overOrUnder = didOverflow ? "over" : "under";
2312 : logger << "BUG: Cache " << overOrUnder << "flow:"
2313 0 : << " nRows=" << msCache.nRows
2314 : << " != nPixels=" << msCache.pixels.size()
2315 0 : << " for: " << msCache.msPath
2316 0 : << " with selection: " << msCache.msTableName
2317 0 : << LogIO::EXCEPTION;
2318 : }
2319 : }
2320 : }
2321 : // Reset state
2322 228 : isOpened = false;
2323 228 : accessMode = AccessMode::READ;
2324 228 : canRead = false;
2325 228 : canWrite = false;
2326 228 : }
2327 :
2328 : Bool
2329 352465 : SDGrid::Cache::isReadable() const {
2330 352465 : return isOpened and canRead;
2331 : }
2332 :
2333 : Bool
2334 176525 : SDGrid::Cache::isWriteable() const {
2335 176525 : return isOpened and canWrite;
2336 : }
2337 :
2338 : Bool
2339 228 : SDGrid::Cache::isEmpty() const {
2340 228 : return msCaches.empty();
2341 : }
2342 :
2343 : void
2344 0 : SDGrid::Cache::clear() {
2345 0 : msCaches.clear();
2346 0 : }
2347 :
2348 129 : SDGrid::Cache::MsCache::MsCache(const String& msPathIn, const String& msTableNameIn, rownr_t nRowsIn)
2349 : : msPath {msPathIn},
2350 : msTableName {msTableNameIn},
2351 129 : nRows {nRowsIn}
2352 : {
2353 129 : pixels.reserve(nRows);
2354 129 : }
2355 :
2356 129 : Bool SDGrid::Cache::MsCache::isConsistent() const {
2357 129 : return pixels.size() == nRows;
2358 : }
2359 :
2360 : void
2361 258 : SDGrid::Cache::newMS(const MeasurementSet& ms) {
2362 516 : LogIO os(LogOrigin(className(),"newMS"));
2363 258 : const auto msPath = ms.getPartNames()[0];
2364 :
2365 258 : if (isWriteable()) {
2366 129 : os << "Will compute and cache spectra pixels coordinates for: " << msPath << LogIO::POST;
2367 129 : msCaches.emplace_back(msPath, ms.tableName(), ms.nrow());
2368 129 : msPixels = &(msCaches.back().pixels);
2369 129 : return;
2370 : }
2371 :
2372 129 : if (isReadable()) {
2373 129 : os << "Will load cached spectra pixels coordinates for: " << msPath << LogIO::POST;
2374 129 : if (msCacheReadIterator == msCaches.cend()) {
2375 0 : os << "BUG! Cached data missing for: " << msPath << LogIO::EXCEPTION;
2376 : }
2377 129 : const auto & pixelsToLoad = msCacheReadIterator->pixels;
2378 129 : if (pixelsToLoad.size() != ms.nrow()) {
2379 : os << "BUG! Cached data size mismatch for: " << msPath
2380 : << " : nRows: " << ms.nrow()
2381 0 : << " nCachedRows: " << pixelsToLoad.size() << LogIO::EXCEPTION;
2382 : }
2383 129 : if ( msCacheReadIterator->msTableName != ms.tableName()) {
2384 : os << "BUG! Cached data was probably for a different selection of: " << msPath
2385 : << " current selection: " << ms.tableName()
2386 0 : << " cached selection: " << msCacheReadIterator->msTableName << LogIO::EXCEPTION;
2387 : }
2388 129 : pixelReadIterator = pixelsToLoad.cbegin();
2389 129 : ++msCacheReadIterator;
2390 : }
2391 : }
2392 :
2393 : void
2394 176039 : SDGrid::Cache::storeRowPixel() {
2395 176039 : msPixels->emplace_back(
2396 176039 : inputPixel.xy[0],
2397 176039 : inputPixel.xy[1],
2398 176039 : inputPixel.isValid
2399 : );
2400 176039 : }
2401 :
2402 : void
2403 176039 : SDGrid::Cache::loadRowPixel() {
2404 176039 : const auto & pixel = *pixelReadIterator;
2405 176039 : outputPixel.xy[0] = pixel.x;
2406 176039 : outputPixel.xy[1] = pixel.y;
2407 176039 : outputPixel.isValid = pixel.isValid;
2408 176039 : ++pixelReadIterator;
2409 176039 : }
2410 :
2411 228 : SDGrid::CacheManager::CacheManager(Cache& cacheIn,
2412 228 : Bool onDutyIn, Cache::AccessMode accessModeIn)
2413 : : cache {cacheIn},
2414 : onDuty {onDutyIn},
2415 228 : accessMode {accessModeIn}
2416 : {
2417 228 : if (onDuty) {
2418 228 : cache.open(accessMode);
2419 : }
2420 228 : }
2421 :
2422 456 : SDGrid::CacheManager::~CacheManager()
2423 : {
2424 228 : if (onDuty) {
2425 228 : cache.close();
2426 : }
2427 228 : }
2428 :
2429 179047 : SDGrid::CacheWriter::CacheWriter(Cache& cacheIn,
2430 179047 : Bool onDutyIn)
2431 : : cache {cacheIn},
2432 179047 : onDuty {onDutyIn}
2433 179047 : {}
2434 :
2435 358094 : SDGrid::CacheWriter::~CacheWriter()
2436 : {
2437 179047 : if (onDuty) {
2438 176039 : cache.storeRowPixel();
2439 : }
2440 179047 : }
2441 :
2442 114 : const String & SDGrid::toString(const ConvertFirst convertFirst) {
2443 : static const std::array<String,3> name {
2444 : "NEVER",
2445 : "ALWAYS",
2446 : "AUTO"
2447 114 : };
2448 114 : switch(convertFirst){
2449 114 : case ConvertFirst::NEVER:
2450 : case ConvertFirst::ALWAYS:
2451 : case ConvertFirst::AUTO:
2452 114 : return name[static_cast<size_t>(convertFirst)];
2453 0 : default:
2454 0 : String errMsg {"Illegal ConvertFirst enum: "};
2455 0 : errMsg += String::toString(static_cast<Int>(convertFirst));
2456 : throw AipsError(
2457 : errMsg,
2458 : __FILE__,
2459 : __LINE__,
2460 : AipsError::Category::INVALID_ARGUMENT
2461 0 : );
2462 : // Avoid potential compiler warning
2463 : return name[static_cast<size_t>(ConvertFirst::NEVER)];
2464 : }
2465 : }
2466 :
2467 114 : SDGrid::ConvertFirst SDGrid::fromString(const String & name) {
2468 : static const std::array<ConvertFirst,3> schemes {
2469 : ConvertFirst::NEVER,
2470 : ConvertFirst::ALWAYS,
2471 : ConvertFirst::AUTO
2472 : };
2473 114 : for ( const auto scheme : schemes ) {
2474 114 : if (name == toString(scheme)) return scheme;
2475 : }
2476 0 : String errMsg {"Illegal ConvertFirst name: "};
2477 0 : errMsg += name;
2478 : throw AipsError(
2479 : errMsg,
2480 : __FILE__,
2481 : __LINE__,
2482 : AipsError::Category::INVALID_ARGUMENT
2483 0 : );
2484 : // Avoid potential compiler warning
2485 : return ConvertFirst::NEVER;
2486 : }
2487 :
2488 114 : void SDGrid::setConvertFirst(const casacore::String &name) {
2489 114 : processingScheme = fromString(name);
2490 114 : }
2491 :
2492 258 : Bool SDGrid::mustConvertPointingColumn(
2493 : const MeasurementSet &ms
2494 : ) {
2495 : const auto haveCachedSpectraPixelCoordinates =
2496 258 : cacheIsEnabled and cache.isReadable();
2497 258 : if (haveCachedSpectraPixelCoordinates) return False;
2498 :
2499 129 : switch(processingScheme){
2500 0 : case ConvertFirst::ALWAYS: return True;
2501 129 : case ConvertFirst::NEVER: return False;
2502 0 : case ConvertFirst::AUTO:
2503 : {
2504 0 : const auto nPointings = ms.pointing().nrow();
2505 0 : const auto nSelectedDataRows = ms.nrow();
2506 0 : return nSelectedDataRows > nPointings ? True : False;
2507 : }
2508 0 : default:
2509 0 : String errMsg {"Unexpected invalid state: "};
2510 0 : errMsg += "ConvertFirst processingScheme=";
2511 0 : errMsg += String::toString<Int>(static_cast<Int>(processingScheme));
2512 0 : errMsg += " ms=" + ms.tableName();
2513 : throw AipsError(
2514 : errMsg,
2515 : __FILE__,
2516 : __LINE__,
2517 : AipsError::Category::GENERAL
2518 0 : );
2519 : }
2520 : // Avoid potential compiler warning
2521 : return False;
2522 : }
2523 :
2524 258 : void SDGrid::handleNewMs(
2525 : ROVisibilityIterator &vi,
2526 : const ImageInterface<Complex>& image) {
2527 :
2528 : // Synchronize spatial coordinates cache
2529 258 : if (cacheIsEnabled) cache.newMS(vi.ms());
2530 :
2531 : // Handle interpolate-convert processing scheme
2532 258 : if (mustConvertPointingColumn(vi.ms())) {
2533 0 : const auto columnEnum = MSPointing::columnType(pointingDirCol_p);
2534 : const auto refType =
2535 0 : image.coordinates().directionCoordinate().directionType();
2536 0 : convertPointingColumn(vi.ms(), columnEnum, refType);
2537 : }
2538 : else {
2539 258 : ramPointingTable = MSPointing();
2540 258 : ramPointingColumnsPtr.reset();
2541 : }
2542 258 : }
2543 :
2544 : namespace convert_pointing_column_helpers {
2545 : using DirectionComputer = std::function<MDirection (rownr_t)>;
2546 :
2547 : DirectionComputer
2548 0 : metaDirectionComputer(
2549 : const MSPointingColumns &pointingColumns,
2550 : MSPointing::PredefinedColumns columnEnum) {
2551 : using std::placeholders::_1;
2552 : using Column = MSPointing::PredefinedColumns;
2553 0 : const auto & encoderDirections = pointingColumns.encoderMeas();
2554 0 : switch(columnEnum) {
2555 0 : case Column::DIRECTION:
2556 0 : return std::bind(
2557 0 : &MSPointingColumns::directionMeas,
2558 0 : &pointingColumns,
2559 : _1,
2560 0 : Double {0.0}
2561 0 : );
2562 : break;
2563 0 : case Column::TARGET:
2564 0 : return std::bind(
2565 0 : &MSPointingColumns::targetMeas,
2566 0 : &pointingColumns,
2567 : _1,
2568 0 : Double {0.0}
2569 0 : );
2570 : break;
2571 0 : case Column::SOURCE_OFFSET:
2572 0 : return std::bind(
2573 0 : &MSPointingColumns::sourceOffsetMeas,
2574 0 : &pointingColumns,
2575 : _1,
2576 0 : Double {0.0}
2577 0 : );
2578 : break;
2579 0 : case Column::POINTING_OFFSET:
2580 0 : return std::bind(
2581 0 : &MSPointingColumns::pointingOffsetMeas,
2582 0 : &pointingColumns,
2583 : _1,
2584 0 : Double {0.0}
2585 0 : );
2586 : break;
2587 0 : case Column::ENCODER:
2588 0 : return std::bind(
2589 0 : &ScalarMeasColumn<MDirection>::operator(),
2590 0 : &encoderDirections,
2591 : _1
2592 0 : );
2593 : break;
2594 0 : default:
2595 : throw AipsError(
2596 0 : String("Illegal Pointing Column Enum: " + String::toString(columnEnum)),
2597 : AipsError::INVALID_ARGUMENT
2598 0 : );
2599 : }
2600 : }
2601 :
2602 : // DirectionArchiver
2603 : class DirectionArchiver {
2604 : public:
2605 : virtual void put(rownr_t row, const MDirection & dir) = 0;
2606 : };
2607 :
2608 : // Derived Templated Class,
2609 : // implementing the shared constructor.
2610 : template<typename ColumnType, typename CellType>
2611 : class DirectionArchiver_ : public DirectionArchiver {
2612 : public:
2613 0 : DirectionArchiver_(ColumnType &columnIn)
2614 0 : : column {columnIn}
2615 0 : {}
2616 : void put(rownr_t row, const MDirection & dir);
2617 : private:
2618 : ColumnType & column;
2619 : };
2620 :
2621 : // Specializations of "put" member
2622 : template<>
2623 0 : void DirectionArchiver_<MDirection::ScalarColumn, MDirection>::put(
2624 : rownr_t row, const MDirection &dir) {
2625 0 : column.put(row, dir);
2626 0 : }
2627 :
2628 : template<>
2629 0 : void DirectionArchiver_<MDirection::ArrayColumn, Array<MDirection>>::put(
2630 : rownr_t row, const MDirection &dir) {
2631 0 : column.put(row, Vector<MDirection> {dir});
2632 0 : }
2633 :
2634 : // Template instantiations for the types we are interested in.
2635 : template class DirectionArchiver_ <MDirection::ArrayColumn, Array<MDirection>>;
2636 : template class DirectionArchiver_ <MDirection::ScalarColumn, MDirection>;
2637 :
2638 : // Aliases
2639 : using ScalarArchiver =
2640 : DirectionArchiver_ <MDirection::ScalarColumn, MDirection>;
2641 : using ArrayArchiver =
2642 : DirectionArchiver_ <MDirection::ArrayColumn, Array<MDirection>>;
2643 :
2644 :
2645 : struct ArchiverFactory {
2646 : static DirectionArchiver * createArchiver(
2647 : TableMeasColumn & column
2648 : );
2649 : };
2650 :
2651 : DirectionArchiver *
2652 0 : ArchiverFactory::createArchiver(TableMeasColumn &column) {
2653 : try {
2654 0 : auto & scalarColumn = dynamic_cast<MDirection::ScalarColumn &>(column);
2655 0 : return new ScalarArchiver(scalarColumn);
2656 : }
2657 0 : catch(std::bad_cast & exception) {
2658 0 : auto & arrayColumn = dynamic_cast<MDirection::ArrayColumn &>(column);
2659 0 : return new ArrayArchiver(arrayColumn);
2660 : }
2661 : }
2662 :
2663 : TableMeasColumn &
2664 0 : columnData(
2665 : MSPointingColumns & pointingColumns,
2666 : MSPointing::PredefinedColumns columnEnum) {
2667 : using Column = MSPointing::PredefinedColumns;
2668 0 : switch(columnEnum){
2669 : // Array Columns
2670 0 : case Column::DIRECTION:
2671 0 : return pointingColumns.directionMeasCol();
2672 0 : case Column::TARGET:
2673 0 : return pointingColumns.targetMeasCol();
2674 0 : case Column::SOURCE_OFFSET:
2675 0 : return pointingColumns.sourceOffsetMeasCol();
2676 0 : case Column::POINTING_OFFSET:
2677 0 : return pointingColumns.pointingOffsetMeasCol();
2678 : // Scalar Column
2679 0 : case Column::ENCODER:
2680 0 : return pointingColumns.encoderMeas();
2681 0 : default:
2682 : {
2683 0 : LogIO logger {LogOrigin {"columnData"} };
2684 : logger << LogIO::EXCEPTION
2685 : << "Expected a column of directions, got: " << MSPointing::columnName(columnEnum)
2686 0 : << LogIO::POST;
2687 :
2688 : // This is just to silence the following compiler warning:
2689 : // warning: control reaches end of non-void function [-Wreturn-type]
2690 0 : return pointingColumns.directionMeasCol();
2691 : }
2692 : }
2693 : }
2694 :
2695 : } // namespace convert_pointing_column_helpers
2696 : using namespace convert_pointing_column_helpers;
2697 :
2698 0 : void SDGrid::convertPointingColumn(
2699 : const MeasurementSet & ms,
2700 : const MSPointingEnums::PredefinedColumns columnEnum,
2701 : const MDirection::Types refType) {
2702 :
2703 0 : LogIO logger {LogOrigin {"SDGrid", "convertPointingColumn"}};
2704 0 : logger << "Start" << LogIO::POST;
2705 :
2706 0 : initRamPointingTable(ms.pointing(), columnEnum, refType);
2707 :
2708 : // Setup helper tools
2709 : // ---- Conversion Tools
2710 0 : MeasFrame pointingMeasurements;
2711 0 : MDirection::Convert convertToImageDirectionRef;
2712 0 : std::tie(
2713 : pointingMeasurements,
2714 : convertToImageDirectionRef
2715 0 : ) = setupConversionTools(ms, refType);
2716 :
2717 : // ---- Direction Computer
2718 0 : MSPointingColumns pointingColumns {ms.pointing()};
2719 : auto userSpecifiedDirection =
2720 0 : metaDirectionComputer(pointingColumns, columnEnum);
2721 :
2722 : // ---- Direction Archiver
2723 0 : MSPointingColumns ramPointingColumns {ramPointingTable};
2724 : std::unique_ptr<DirectionArchiver> correspondingRamPointingColumn {
2725 : ArchiverFactory::createArchiver(
2726 : columnData(
2727 : ramPointingColumns,
2728 : columnEnum
2729 : )
2730 : )
2731 0 : };
2732 :
2733 : // Convert directions stored in user-specified
2734 : // POINTING column (of some direction),
2735 : // and store the result into the corresponding column
2736 : // of the RAM POINTING table
2737 0 : const auto & epoch = pointingColumns.timeMeas();
2738 0 : const auto & antennaId = pointingColumns.antennaId();
2739 :
2740 0 : MSAntennaColumns antennaColumns(ms.antenna());
2741 0 : const auto & antennaPosition = antennaColumns.positionMeas();
2742 :
2743 : // Main loop control
2744 0 : const auto pointingRows = ms.pointing().nrow();
2745 0 : constexpr Int invalidAntennaId = -1;
2746 0 : auto previousPointing_AntennaId {invalidAntennaId};
2747 :
2748 : // Main loop
2749 0 : for (rownr_t pointingRow = 0; pointingRow < pointingRows ; ++pointingRow){
2750 : // Collect pointing measurements
2751 0 : const auto pointing_Epoch = epoch(pointingRow);
2752 0 : pointingMeasurements.resetEpoch(pointing_Epoch);
2753 :
2754 0 : const auto pointing_AntennaId = antennaId(pointingRow);
2755 0 : const auto antennaChanged =
2756 : ( pointing_AntennaId != previousPointing_AntennaId );
2757 0 : if (antennaChanged) {
2758 : const auto new_AntennaPosition =
2759 0 : antennaPosition(pointing_AntennaId);
2760 0 : pointingMeasurements.resetPosition(new_AntennaPosition);
2761 0 : previousPointing_AntennaId = pointing_AntennaId;
2762 : }
2763 :
2764 : // Convert
2765 : const auto convertedDirection =
2766 : convertToImageDirectionRef(
2767 0 : userSpecifiedDirection(pointingRow)
2768 0 : );
2769 :
2770 : // Store
2771 0 : correspondingRamPointingColumn->put(
2772 : pointingRow,
2773 : convertedDirection
2774 0 : );
2775 : }
2776 :
2777 0 : logger << "Done" << LogIO::POST;
2778 0 : }
2779 :
2780 0 : void SDGrid::initRamPointingTable(
2781 : const MSPointing & pointingTable,
2782 : const MSPointingEnums::PredefinedColumns columnEnum,
2783 : const MDirection::Types refType) {
2784 :
2785 0 : LogIO logger {LogOrigin {"SDGrid", "initRamPointingTable"}};
2786 0 : logger << "Start" << LogIO::POST;
2787 :
2788 0 : constexpr auto doNotCopyRows = True;
2789 0 : ramPointingTable = pointingTable.copyToMemoryTable(
2790 0 : pointingTable.tableName() +
2791 0 : "." + MSPointing::columnName(columnEnum) +
2792 0 : "." + MDirection::showType(refType),
2793 : doNotCopyRows
2794 0 : );
2795 0 : ramPointingColumnsPtr.reset( new MSPointingColumns {ramPointingTable});
2796 :
2797 0 : MSPointingColumns pointingColumns {pointingTable};
2798 :
2799 0 : ramPointingColumnsPtr->setDirectionRef(refType);
2800 0 : ramPointingColumnsPtr->setEncoderDirectionRef(refType);
2801 :
2802 0 : ramPointingTable.addRow(pointingTable.nrow());
2803 :
2804 0 : ramPointingColumnsPtr->antennaId().putColumn(pointingColumns.antennaId());
2805 0 : ramPointingColumnsPtr->time().putColumn(pointingColumns.time());
2806 0 : ramPointingColumnsPtr->interval().putColumn(pointingColumns.interval());
2807 0 : ramPointingColumnsPtr->numPoly().fillColumn(0);
2808 :
2809 0 : logger << "Done" << LogIO::POST;
2810 0 : }
2811 :
2812 : std::pair<MeasFrame,MDirection::Convert>
2813 0 : SDGrid::setupConversionTools(
2814 : const MeasurementSet & ms,
2815 : const casacore::MDirection::Types refType) {
2816 :
2817 0 : LogIO logger {LogOrigin {"SDGrid", "setupConversionTools"}};
2818 0 : logger << "Start" << LogIO::POST;
2819 :
2820 0 : MSPointingColumns pointingColumns(ms.pointing());
2821 0 : MSAntennaColumns antennaColumns(ms.antenna());
2822 :
2823 0 : auto firstPointing_Epoch = pointingColumns.timeMeas()(0);
2824 0 : auto firstPointing_AntennaId = pointingColumns.antennaId()(0);
2825 0 : auto firstPointing_AntennaPosition = antennaColumns.positionMeas()(
2826 : static_cast<rownr_t>(firstPointing_AntennaId)
2827 0 : );
2828 :
2829 0 : MeasFrame measFrame(firstPointing_Epoch, firstPointing_AntennaPosition);
2830 :
2831 : // Direction Conversion Machine
2832 0 : MDirection::Ref dstRef(refType, measFrame);
2833 0 : auto firstPointing_Direction = pointingColumns.directionMeas(0);
2834 0 : MDirection::Convert convert(firstPointing_Direction, dstRef);
2835 :
2836 : // ---- Perform 1 dummy conversion at a time far before
2837 : // the first pointing time, so that when we will next convert
2838 : // the first user-specified direction,
2839 : // we are sure that values cached in static variables
2840 : // of casacore functions like dUT1() will be cleared
2841 : static const MVEpoch oneYear {
2842 0 : Quantity {
2843 : 365,
2844 0 : Unit { "d" }
2845 : }
2846 0 : };
2847 : const MEpoch dummy_Epoch {
2848 0 : firstPointing_Epoch.getValue() - oneYear,
2849 0 : firstPointing_Epoch.getRef()
2850 0 : };
2851 0 : measFrame.resetEpoch(dummy_Epoch);
2852 0 : const auto dummy_Direction = convert();
2853 :
2854 0 : logger << "Done" << LogIO::POST;
2855 0 : return std::make_pair(measFrame, convert);
2856 : }
2857 :
2858 : } // casa namespace
|
