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 0 : ChronoStat::ChronoStat(const string & name)
93 : : name_ {name},
94 : started_ {false},
95 : n_laps_ {0},
96 : n_overflows_ {0},
97 : n_underflows_ {0},
98 0 : laps_sum_ {Duration::zero()},
99 0 : laps_min_ {Duration::max()},
100 0 : laps_max_ {Duration::min()}
101 0 : {}
102 :
103 0 : const std::string& ChronoStat::name() const {
104 0 : return name_;
105 : }
106 :
107 0 : void ChronoStat::setName(const std::string& name) {
108 0 : name_ = name;
109 0 : }
110 :
111 0 : void ChronoStat::start() {
112 0 : if (not started_) {
113 0 : started_ = true;
114 0 : ++n_laps_;
115 0 : lap_start_time_ = Clock::now();
116 : } else {
117 0 : ++n_overflows_;
118 : }
119 0 : }
120 0 : void ChronoStat::stop() {
121 0 : if (started_) {
122 0 : auto lap_duration = Clock::now() - lap_start_time_;
123 0 : laps_sum_ += lap_duration;
124 0 : started_ = false;
125 0 : if (lap_duration < laps_min_) laps_min_ = lap_duration;
126 0 : if (lap_duration > laps_max_) laps_max_ = lap_duration;
127 : } else {
128 0 : ++n_underflows_;
129 : }
130 0 : }
131 :
132 0 : ChronoStat::Duration ChronoStat::lapsSum() const {
133 0 : return laps_sum_;
134 : }
135 :
136 0 : ChronoStat::Duration ChronoStat::lapsMin() const {
137 0 : return n_laps_ > 0 ? laps_min_ : Duration::zero();
138 : }
139 :
140 0 : ChronoStat::Duration ChronoStat::lapsMax() const {
141 0 : return n_laps_ > 0 ? laps_max_ : Duration::zero();
142 : }
143 :
144 0 : ChronoStat::Duration ChronoStat::lapsMean() const {
145 0 : return n_laps_ > 0 ? laps_sum_ / n_laps_ : Duration::zero();
146 : }
147 :
148 0 : unsigned int ChronoStat::lapsCount() const {
149 0 : return n_laps_;
150 : }
151 :
152 0 : bool ChronoStat::isEmpty() const {
153 0 : return n_laps_ == 0;
154 : }
155 :
156 0 : unsigned int ChronoStat::nOverflows() const {
157 0 : return n_overflows_;
158 : }
159 :
160 0 : unsigned int ChronoStat::nUnderflows() const {
161 0 : return n_underflows_;
162 : }
163 :
164 0 : std::string ChronoStat::quote(const std::string& s) const {
165 0 : return "\"" + s + "\"";
166 : }
167 :
168 0 : std::string ChronoStat::json() const {
169 0 : std::ostringstream os;
170 0 : os << quote(name()) << ": {"
171 0 : << quote("sum") << ": " << lapsSum().count()
172 0 : << " ," << quote("count") << ": " << lapsCount()
173 0 : << " ," << quote("min") << ": " << lapsMin().count()
174 0 : << " ," << quote("mean") << ": " << lapsMean().count()
175 0 : << " ," << quote("max") << ": " << lapsMax().count()
176 0 : << " ," << quote("overflows") << ": " << nOverflows()
177 0 : << " ," << quote("underflows") << ": " << nUnderflows()
178 0 : << "}";
179 0 : 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 0 : StartStop::StartStop(ChronoStat &c)
196 0 : : c_ {c} {
197 0 : c_.start();
198 0 : }
199 :
200 0 : StartStop::~StartStop() {
201 0 : c_.stop();
202 0 : }
203 :
204 : } // namespace sdgrid_perfs
205 :
206 : using namespace sdgrid_perfs;
207 :
208 : #endif
209 :
210 0 : SDGrid::SDGrid(SkyJones& sj, Int icachesize, Int itilesize,
211 0 : 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 0 : 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 0 : cacheIsEnabled {false}
221 : {
222 0 : lastIndex_p=0;
223 0 : initPerfs();
224 0 : }
225 :
226 0 : SDGrid::SDGrid(MPosition& mLocation, SkyJones& sj, Int icachesize, Int itilesize,
227 0 : 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 0 : 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 0 : cacheIsEnabled {false}
237 : {
238 0 : mLocation_p=mLocation;
239 0 : lastIndex_p=0;
240 0 : initPerfs();
241 0 : }
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 0 : SDGrid::SDGrid(MPosition &mLocation, Int icachesize, Int itilesize,
260 0 : 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 0 : 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 0 : cacheIsEnabled {false}
271 : {
272 0 : mLocation_p=mLocation;
273 0 : lastIndex_p=0;
274 0 : initPerfs();
275 0 : }
276 :
277 0 : SDGrid::SDGrid(MPosition &mLocation, Int icachesize, Int itilesize,
278 : String iconvType, Float truncate, Float gwidth, Float jwidth,
279 0 : 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 0 : 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 0 : cacheIsEnabled {false}
289 : {
290 0 : mLocation_p=mLocation;
291 0 : lastIndex_p=0;
292 0 : initPerfs();
293 0 : }
294 :
295 : //----------------------------------------------------------------------
296 :
297 0 : void SDGrid::initPerfs() {
298 : #if defined(SDGRID_PERFS)
299 0 : cNextChunk.setName("iterateNextChunk");
300 0 : cMatchAllSpwChans.setName("matchAllSpwChans");
301 0 : cMatchChannel.setName("matchChannel");
302 0 : cPickWeights.setName("pickWeights");
303 0 : cInterpolateFrequencyToGrid.setName("interpolateFrequencyToGrid");
304 0 : cSearchValidPointing.setName("searchValidPointing");
305 0 : cComputeSplines.setName("computeSplines");
306 0 : cResetFrame.setName("resetFrame");
307 0 : cInterpolateDirection.setName("interpolateDirection");
308 0 : cConvertDirection.setName("convertDirection");
309 0 : cComputeDirectionPixel.setName("computeDirectionPixel");
310 0 : cHandleMovingSource.setName("handleMovingSource");
311 0 : cGridData.setName("gridData");
312 : #endif
313 0 : }
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 0 : String SDGrid::name() const{
364 0 : return String("SDGrid");
365 : }
366 :
367 : void
368 0 : SDGrid::setEnableCache(Bool doEnable) {
369 0 : cacheIsEnabled = doEnable;
370 0 : }
371 :
372 : //----------------------------------------------------------------------
373 : // Odds are that it changed.....
374 0 : Bool SDGrid::changed(const VisBuffer& /*vb*/) {
375 0 : 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 0 : void SDGrid::init() {
412 :
413 0 : logIO() << LogOrigin("SDGrid", "init") << LogIO::NORMAL;
414 :
415 0 : ok();
416 :
417 0 : isTiled = false;
418 0 : nx = image->shape()(0);
419 0 : ny = image->shape()(1);
420 0 : npol = image->shape()(2);
421 0 : nchan = image->shape()(3);
422 :
423 0 : 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 0 : if (imageCache) delete imageCache;
430 0 : imageCache = nullptr;
431 :
432 0 : convType = downcase(convType);
433 0 : logIO() << "Convolution function : " << convType << LogIO::DEBUG1 << LogIO::POST;
434 :
435 : // Compute convolution function
436 0 : if (convType=="pb") {
437 : // Do nothing
438 : }
439 0 : else if (convType=="box") {
440 0 : convSupport=(userSetSupport_p >= 0) ? userSetSupport_p : 0;
441 0 : logIO() << "Support : " << convSupport << " pixels" << LogIO::POST;
442 0 : convSampling=100;
443 0 : convSize=convSampling*(2*convSupport+2);
444 0 : convFunc.resize(convSize);
445 0 : convFunc=0.0;
446 0 : for (Int i=0;i<convSize/2;i++) {
447 0 : convFunc(i)=1.0;
448 : }
449 : }
450 0 : else if (convType=="sf") {
451 0 : convSupport=(userSetSupport_p >= 0) ? userSetSupport_p : 3;
452 0 : logIO() << "Support : " << convSupport << " pixels" << LogIO::POST;
453 0 : convSampling=100;
454 0 : convSize=convSampling*(2*convSupport+2);
455 0 : convFunc.resize(convSize);
456 0 : convFunc=0.0;
457 0 : for (Int i=0;i<convSampling*convSupport;i++) {
458 0 : Double nu=Double(i)/Double(convSupport*convSampling);
459 : Double val;
460 0 : grdsf(&nu, &val);
461 0 : convFunc(i)=(1.0-nu*nu)*val;
462 : }
463 : }
464 0 : else if (convType=="gauss") {
465 : // default is b=1.0 (Mangum et al. 2007)
466 0 : Double hwhm=(gwidth_p > 0.0) ? Double(gwidth_p) : sqrt(log(2.0));
467 0 : Float truncate=(truncate_p >= 0.0) ? truncate_p : 3.0*hwhm;
468 0 : convSampling=100;
469 0 : Int itruncate=(Int)(truncate*Double(convSampling)+0.5);
470 0 : logIO() << LogIO::DEBUG1 << "hwhm=" << hwhm << LogIO::POST;
471 0 : logIO() << LogIO::DEBUG1 << "itruncate=" << itruncate << LogIO::POST;
472 : //convSupport=(Int)(truncate+0.5);
473 0 : convSupport = (Int)(truncate);
474 0 : convSupport += (((truncate-(Float)convSupport) > 0.0) ? 1 : 0);
475 0 : convSize=convSampling*(2*convSupport+2);
476 0 : convFunc.resize(convSize);
477 0 : convFunc=0.0;
478 : Double val, x;
479 0 : for (Int i = 0 ; i <= itruncate ; i++) {
480 0 : x = Double(i)/Double(convSampling);
481 0 : grdgauss(&hwhm, &x, &val);
482 0 : convFunc(i)=val;
483 : }
484 : }
485 0 : else if (convType=="gjinc") {
486 : // default is b=2.52, c=1.55 (Mangum et al. 2007)
487 0 : Double hwhm = (gwidth_p > 0.0) ? Double(gwidth_p) : sqrt(log(2.0))*2.52;
488 0 : Double c = (jwidth_p > 0.0) ? Double(jwidth_p) : 1.55;
489 : //Float truncate = truncate_p;
490 0 : convSampling = 100;
491 0 : Int itruncate=(Int)(truncate_p*Double(convSampling)+0.5);
492 0 : logIO() << LogIO::DEBUG1 << "hwhm=" << hwhm << LogIO::POST;
493 0 : logIO() << LogIO::DEBUG1 << "c=" << c << LogIO::POST;
494 0 : logIO() << LogIO::DEBUG1 << "itruncate=" << itruncate << LogIO::POST;
495 : //convSupport=(truncate_p >= 0.0) ? (Int)(truncate_p+0.5) : (Int)(2*c+0.5);
496 0 : Float convSupportF = (truncate_p >= 0.0) ? truncate_p : (2*c);
497 0 : convSupport = (Int)convSupportF;
498 0 : convSupport += (((convSupportF-(Float)convSupport) > 0.0) ? 1 : 0);
499 0 : convSize=convSampling*(2*convSupport+2);
500 0 : convFunc.resize(convSize);
501 0 : convFunc=0.0;
502 : //UNUSED: Double r;
503 : Double x, val1, val2;
504 0 : Int normalize = 1;
505 0 : 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 0 : for (Int i=0;i<convSize;i++) {
517 0 : x = Double(i) / Double(convSampling);
518 : //r = Double(i) / (Double(hwhm)*Double(convSampling));
519 0 : grdjinc1(&c, &x, &normalize, &val2);
520 0 : if (val2 <= 0.0) {
521 : logIO() << LogIO::DEBUG1
522 : << "convFunc is automatically truncated at radius "
523 0 : << x << LogIO::POST;
524 0 : break;
525 : }
526 0 : grdgauss(&hwhm, &x, &val1);
527 0 : 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 0 : if (wImage) delete wImage;
540 0 : wImage = new TempImage<Float>(image->shape(), image->coordinates());
541 :
542 0 : }
543 :
544 0 : void SDGrid::collectPerfs(){
545 : #if defined(SDGRID_PERFS)
546 0 : 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 0 : };
562 : os << "PERFS<SDGRID> "
563 : << "{ \"note\": \"sum, min, mean, max are in units of nanoseconds.\""
564 : << ", \"probes\": "
565 0 : << "{ " << join(probes_json.data(), probes_json.size(), ", " ) << " }"
566 : << " }"
567 0 : << LogIO::POST;
568 : #endif
569 0 : }
570 :
571 :
572 : // This is nasty, we should use CountedPointers here.
573 0 : SDGrid::~SDGrid() {
574 : //fclose(pfile);
575 0 : if (imageCache) delete imageCache; imageCache = 0;
576 0 : if (arrayLattice) delete arrayLattice; arrayLattice = 0;
577 0 : if (wImage) delete wImage; wImage = 0;
578 0 : if (wImageCache) delete wImageCache; wImageCache = 0;
579 0 : if (wArrayLattice) delete wArrayLattice; wArrayLattice = 0;
580 0 : if (interpolator) delete interpolator; interpolator = 0;
581 :
582 0 : collectPerfs();
583 0 : }
584 :
585 0 : 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 0 : 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 0 : convSupport=max(128, sj_p->support(vb, coords));
595 :
596 0 : convSampling=100;
597 0 : 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 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
603 0 : AlwaysAssert(directionIndex>=0, AipsError);
604 0 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
605 0 : Vector<Double> sampling;
606 0 : sampling = dc.increment();
607 0 : sampling*=1.0/Double(convSampling);
608 0 : 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 0 : uInt row=0;
614 :
615 : // reset lastAntID_p to use correct antenna position
616 0 : lastAntID_p = -1;
617 :
618 0 : const MSPointingColumns& act_mspc = vb.msColumns().pointing();
619 0 : Bool nullPointingTable=(act_mspc.nrow() < 1);
620 : // uInt pointIndex=getIndex(*mspc, vb.time()(row), vb.timeInterval()(row), vb.antenna1()(row));
621 0 : Int pointIndex=-1;
622 0 : if(!nullPointingTable){
623 : //if(vb.newMS()) This thing is buggy...using msId change
624 0 : if(vb.msId() != msId_p){
625 0 : lastIndex_p=0;
626 0 : if(lastIndexPerAnt_p.nelements() < (size_t)vb.numberAnt()) {
627 0 : lastIndexPerAnt_p.resize(vb.numberAnt());
628 : }
629 0 : lastIndexPerAnt_p=0;
630 0 : msId_p=vb.msId();
631 : }
632 0 : pointIndex=getIndex(act_mspc, vb.time()(row), -1.0, vb.antenna1()(row));
633 0 : if(pointIndex < 0)
634 0 : pointIndex=getIndex(act_mspc, vb.time()(row), vb.timeInterval()(row), vb.antenna1()(row));
635 :
636 : }
637 0 : 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 0 : MEpoch epoch(Quantity(vb.time()(row), "s"));
644 0 : if(!pointingToImage) {
645 0 : lastAntID_p=vb.antenna1()(row);
646 0 : MPosition pos;
647 0 : pos=vb.msColumns().antenna().positionMeas()(lastAntID_p);
648 0 : mFrame_p=MeasFrame(epoch, pos);
649 0 : if(!nullPointingTable){
650 0 : 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 0 : MDirection::Ref outRef(dc.directionType(), mFrame_p);
658 0 : pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
659 :
660 0 : 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 0 : if(!nullPointingTable){
678 0 : worldPosMeas=(*pointingToImage)(directionMeas(act_mspc,pointIndex));
679 : }
680 : else{
681 0 : worldPosMeas=(*pointingToImage)(vb.direction1()(row));
682 : }
683 0 : delete pointingToImage;
684 0 : pointingToImage=0;
685 : }
686 :
687 0 : directionCoord=coords.directionCoordinate(directionIndex);
688 : //make sure we use the same units
689 0 : 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 0 : Vector<Double> const originalReferencePixel = dc.referencePixel();
697 0 : dc.setReferenceValue(worldPosMeas.getAngle().getValue());
698 : //Vector<Double> unitVec(2);
699 : //unitVec=0.0;
700 : //dc.setReferencePixel(unitVec);
701 0 : Vector<Double> updatedReferencePixel = dc.referencePixel() - originalReferencePixel;
702 0 : dc.setReferencePixel(updatedReferencePixel);
703 :
704 0 : coords.replaceCoordinate(dc, directionIndex);
705 :
706 0 : IPosition pbShape(4, convSize, 2, 1, 1);
707 0 : IPosition start(4, 0, 0, 0, 0);
708 :
709 0 : TempImage<Complex> onedPB(pbShape, coords);
710 :
711 0 : onedPB.set(Complex(1.0, 0.0));
712 :
713 0 : AlwaysAssert(sj_p, AipsError);
714 0 : sj_p->apply(onedPB, onedPB, vb, 0);
715 :
716 0 : IPosition pbSlice(4, convSize, 1, 1, 1);
717 0 : Vector<Float> tempConvFunc=real(onedPB.getSlice(start, pbSlice, true));
718 : // Find number of significant points
719 0 : uInt cfLen=0;
720 0 : for(uInt i=0;i<tempConvFunc.nelements();++i) {
721 0 : if(tempConvFunc(i)<1e-3) break;
722 0 : ++cfLen;
723 : }
724 0 : 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 0 : Vector<Float> trimConvFunc=tempConvFunc(Slice(0,cfLen-1,1));
731 :
732 : // Now fill in the convolution function vector
733 0 : convSupport=cfLen/convSampling;
734 0 : convSize=convSampling*(2*convSupport+2);
735 0 : convFunc.resize(convSize);
736 0 : convFunc=0.0;
737 0 : convFunc(Slice(0,cfLen-1,1))=trimConvFunc(Slice(0,cfLen-1,1));
738 :
739 :
740 0 : }
741 :
742 : // Initialize for a transform from the Sky domain. This means that
743 : // we grid-correct, and FFT the image
744 0 : void SDGrid::initializeToVis(ImageInterface<Complex>& iimage,
745 : const VisBuffer& vb)
746 : {
747 0 : image=&iimage;
748 :
749 0 : ok();
750 :
751 0 : init();
752 :
753 0 : if(convType=="pb") {
754 0 : 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 0 : msId_p = -1;
760 0 : lastAntID_p = -1;
761 :
762 : // Initialize the maps for polarization and channel. These maps
763 : // translate visibility indices into image indices
764 0 : initMaps(vb);
765 :
766 : // First get the CoordinateSystem for the image and then find
767 : // the DirectionCoordinate
768 0 : CoordinateSystem coords=image->coordinates();
769 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
770 0 : AlwaysAssert(directionIndex>=0, AipsError);
771 0 : directionCoord=coords.directionCoordinate(directionIndex);
772 : /*if((image->shape().product())>cachesize) {
773 : isTiled=true;
774 : }
775 : else {
776 : isTiled=false;
777 : }*/
778 0 : isTiled=false;
779 0 : nx = image->shape()(0);
780 0 : ny = image->shape()(1);
781 0 : npol = image->shape()(2);
782 0 : 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 0 : IPosition gridShape(4, nx, ny, npol, nchan);
794 0 : griddedData.resize(gridShape);
795 0 : griddedData = Complex(0.0);
796 :
797 0 : wGriddedData.resize(gridShape);
798 0 : wGriddedData = 0.0;
799 :
800 0 : if(arrayLattice) delete arrayLattice; arrayLattice=0;
801 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
802 :
803 0 : if(wArrayLattice) delete wArrayLattice; wArrayLattice=0;
804 0 : wArrayLattice = new ArrayLattice<Float>(wGriddedData);
805 0 : wArrayLattice->set(0.0);
806 0 : wLattice=wArrayLattice;
807 :
808 : // Now find the SubLattice corresponding to the image
809 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
810 0 : IPosition stride(4, 1);
811 0 : IPosition trc(blc+image->shape()-stride);
812 0 : LCBox gridBox(blc, trc, gridShape);
813 0 : SubLattice<Complex> gridSub(*arrayLattice, gridBox, true);
814 :
815 : // Do the copy
816 0 : gridSub.copyData(*image);
817 :
818 0 : lattice=arrayLattice;
819 : }
820 :
821 0 : AlwaysAssert(lattice, AipsError);
822 0 : AlwaysAssert(wLattice, AipsError);
823 :
824 0 : }
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 0 : void SDGrid::initializeToSky(ImageInterface<Complex>& iimage,
846 : Matrix<Float>& weight, const VisBuffer& vb)
847 : {
848 : // image always points to the image
849 0 : image=&iimage;
850 :
851 0 : ok();
852 :
853 0 : init();
854 :
855 0 : 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 0 : msId_p = -1;
860 0 : lastAntID_p = -1;
861 :
862 : // Initialize the maps for polarization and channel.
863 : // These maps translate visibility indices into image indices
864 0 : initMaps(vb);
865 :
866 : // No longer using isTiled
867 0 : isTiled=false;
868 0 : nx = image->shape()(0);
869 0 : ny = image->shape()(1);
870 0 : npol = image->shape()(2);
871 0 : nchan = image->shape()(3);
872 :
873 0 : sumWeight = 0.0;
874 0 : weight.resize(sumWeight.shape());
875 0 : weight = 0.0;
876 :
877 : // First get the CoordinateSystem for the image
878 : // and then find the DirectionCoordinate
879 0 : CoordinateSystem coords = image->coordinates();
880 0 : Int directionIndex = coords.findCoordinate(Coordinate::DIRECTION);
881 0 : AlwaysAssert(directionIndex >= 0, AipsError);
882 0 : directionCoord = coords.directionCoordinate(directionIndex);
883 :
884 : // Initialize for in memory gridding.
885 : // lattice will point to the ArrayLattice
886 0 : IPosition gridShape(4, nx, ny, npol, nchan);
887 0 : logIO() << LogOrigin("SDGrid", "initializeToSky", WHERE) << LogIO::DEBUGGING
888 0 : << "gridShape = " << gridShape << LogIO::POST;
889 :
890 0 : griddedData.resize(gridShape);
891 0 : griddedData = Complex(0.0);
892 0 : if (arrayLattice) delete arrayLattice;
893 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
894 0 : lattice = arrayLattice;
895 0 : AlwaysAssert(lattice, AipsError);
896 :
897 0 : wGriddedData.resize(gridShape);
898 0 : wGriddedData=0.0;
899 0 : if (wArrayLattice) delete wArrayLattice;
900 0 : wArrayLattice = new ArrayLattice<Float>(wGriddedData);
901 0 : wLattice = wArrayLattice;
902 0 : AlwaysAssert(wLattice, AipsError);
903 :
904 : // clipping related stuff
905 0 : if (clipminmax_) {
906 0 : gmin_.resize(gridShape);
907 0 : gmin_ = Complex(FLT_MAX);
908 0 : gmax_.resize(gridShape);
909 0 : gmax_ = Complex(-FLT_MAX);
910 0 : wmin_.resize(gridShape);
911 0 : wmin_ = 0.0f;
912 0 : wmax_.resize(gridShape);
913 0 : wmax_ = 0.0f;
914 0 : npoints_.resize(gridShape);
915 0 : npoints_ = 0;
916 : }
917 :
918 : // debug messages
919 0 : LogOrigin msgOrigin("SDGrid", "initializeToSky", WHERE);
920 0 : auto & logger = logIO();
921 0 : logger << msgOrigin << LogIO::DEBUGGING;
922 0 : logger.output() << "clipminmax_ = " << std::boolalpha << clipminmax_
923 0 : << " (" << std::noboolalpha << clipminmax_ << ")";
924 0 : logger << LogIO::POST;
925 0 : if (clipminmax_) {
926 : logger << msgOrigin << LogIO::DEBUGGING
927 : << "will use clipping-capable Fortran gridder ggridsd2 for imaging"
928 0 : << LogIO::POST;
929 : }
930 0 : }
931 :
932 0 : void SDGrid::finalizeToSky()
933 : {
934 0 : if (pointingToImage) delete pointingToImage;
935 0 : pointingToImage = nullptr;
936 0 : }
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 0 : void SDGrid::put(const VisBuffer& vb, Int row, Bool dopsf,
1035 : FTMachine::Type type)
1036 : {
1037 0 : LogIO os(LogOrigin("SDGrid", "put"));
1038 :
1039 0 : gridOk(convSupport);
1040 :
1041 : // Check if ms has changed then cache new spw and chan selection
1042 0 : if (vb.newMS()) {
1043 : #if defined(SDGRID_PERFS)
1044 0 : StartStop trigger(cMatchAllSpwChans);
1045 : #endif
1046 0 : matchAllSpwChans(vb);
1047 0 : lastIndex_p = 0;
1048 0 : if (lastIndexPerAnt_p.nelements() < (size_t)vb.numberAnt()) {
1049 0 : lastIndexPerAnt_p.resize(vb.numberAnt());
1050 : }
1051 0 : 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 0 : if (doConversion_p[vb.spectralWindow()]) {
1057 : #if defined(SDGRID_PERFS)
1058 0 : StartStop trigger(cMatchChannel);
1059 : #endif
1060 0 : matchChannel(vb.spectralWindow(), vb);
1061 : }
1062 : else {
1063 0 : chanMap.resize();
1064 0 : chanMap = multiChanMap_p[vb.spectralWindow()];
1065 : }
1066 :
1067 : // No point in reading data if its not matching in frequency
1068 0 : if (max(chanMap)==-1)
1069 0 : return;
1070 :
1071 0 : Matrix<Float> imagingweight;
1072 0 : pickWeights(vb, imagingweight);
1073 :
1074 0 : if (type==FTMachine::PSF || type==FTMachine::COVERAGE)
1075 0 : dopsf = true;
1076 0 : if (dopsf) type=FTMachine::PSF;
1077 :
1078 0 : Cube<Complex> data;
1079 : //Fortran gridder need the flag as ints
1080 0 : Cube<Int> flags;
1081 0 : Matrix<Float> elWeight;
1082 : {
1083 : #if defined(SDGRID_PERFS)
1084 0 : StartStop trigger(cInterpolateFrequencyToGrid);
1085 : #endif
1086 0 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
1087 : }
1088 :
1089 : Bool iswgtCopy;
1090 : const Float *wgtStorage;
1091 0 : wgtStorage=elWeight.getStorage(iswgtCopy);
1092 : Bool isCopy;
1093 0 : const Complex *datStorage = nullptr;
1094 0 : if (!dopsf)
1095 0 : datStorage = data.getStorage(isCopy);
1096 :
1097 : // If row is -1 then we pass through all rows
1098 : Int startRow, endRow, nRow;
1099 0 : if (row == -1) {
1100 0 : nRow = vb.nRow();
1101 0 : startRow = 0;
1102 0 : endRow = nRow - 1;
1103 : } else {
1104 0 : nRow = 1;
1105 0 : startRow = endRow = row;
1106 : }
1107 :
1108 0 : Vector<Int> rowFlags(vb.flagRow().nelements(),0);
1109 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1110 0 : if (vb.flagRow()(rownr)) rowFlags(rownr) = 1;
1111 : }
1112 :
1113 : // Take care of translation of Bools to Integer
1114 0 : Int idopsf = dopsf ? 1 : 0;
1115 :
1116 : /*if(isTiled) {
1117 : }
1118 : else*/
1119 : {
1120 0 : Matrix<Double> xyPositions(2, endRow - startRow + 1);
1121 0 : xyPositions = -1e9; // make sure failed getXYPos does not fall on grid
1122 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1123 0 : if (getXYPos(vb, rownr)) {
1124 0 : xyPositions(0, rownr) = xyPos(0);
1125 0 : xyPositions(1, rownr) = xyPos(1);
1126 : }
1127 : }
1128 : {
1129 : #if defined(SDGRID_PERFS)
1130 0 : StartStop trigger(cGridData);
1131 : #endif
1132 : Bool del;
1133 : // IPosition s(data.shape());
1134 0 : const IPosition& fs=flags.shape();
1135 0 : std::vector<Int> s(fs.begin(), fs.end());
1136 : Bool datCopy, wgtCopy;
1137 0 : Complex * datStor=griddedData.getStorage(datCopy);
1138 0 : Float * wgtStor=wGriddedData.getStorage(wgtCopy);
1139 :
1140 0 : Bool call_ggridsd = !clipminmax_ || dopsf;
1141 :
1142 0 : if (call_ggridsd) {
1143 :
1144 0 : ggridsd(xyPositions.getStorage(del),
1145 : datStorage,
1146 0 : &s[0],
1147 0 : &s[1],
1148 : &idopsf,
1149 0 : flags.getStorage(del),
1150 0 : rowFlags.getStorage(del),
1151 : wgtStorage,
1152 0 : &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 0 : Complex *gminStor = gmin_.getStorage(gminCopy);
1171 : Bool gmaxCopy;
1172 0 : Complex *gmaxStor = gmax_.getStorage(gmaxCopy);
1173 : Bool wminCopy;
1174 0 : Float *wminStor = wmin_.getStorage(wminCopy);
1175 : Bool wmaxCopy;
1176 0 : Float *wmaxStor = wmax_.getStorage(wmaxCopy);
1177 : Bool npCopy;
1178 0 : Int *npStor = npoints_.getStorage(npCopy);
1179 :
1180 0 : ggridsdclip(xyPositions.getStorage(del),
1181 : datStorage,
1182 0 : &s[0],
1183 0 : &s[1],
1184 : &idopsf,
1185 0 : flags.getStorage(del),
1186 0 : rowFlags.getStorage(del),
1187 : wgtStorage,
1188 0 : &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 0 : gmin_.putStorage(gminStor, gminCopy);
1209 0 : gmax_.putStorage(gmaxStor, gmaxCopy);
1210 0 : wmin_.putStorage(wminStor, wminCopy);
1211 0 : wmax_.putStorage(wmaxStor, wmaxCopy);
1212 0 : npoints_.putStorage(npStor, npCopy);
1213 : }
1214 0 : griddedData.putStorage(datStor, datCopy);
1215 0 : wGriddedData.putStorage(wgtStor, wgtCopy);
1216 : }
1217 : }
1218 0 : if (!dopsf)
1219 0 : data.freeStorage(datStorage, isCopy);
1220 0 : elWeight.freeStorage(wgtStorage,iswgtCopy);
1221 : }
1222 :
1223 0 : void SDGrid::get(VisBuffer& vb, Int row)
1224 : {
1225 0 : LogIO os(LogOrigin("SDGrid", "get"));
1226 :
1227 0 : gridOk(convSupport);
1228 :
1229 : // If row is -1 then we pass through all rows
1230 : Int startRow, endRow, nRow;
1231 0 : if (row==-1) {
1232 0 : nRow=vb.nRow();
1233 0 : startRow=0;
1234 0 : endRow=nRow-1;
1235 0 : 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 0 : 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 0 : if(doConversion_p[vb.spectralWindow()]){
1258 0 : matchChannel(vb.spectralWindow(), vb);
1259 : }
1260 : else{
1261 0 : chanMap.resize();
1262 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
1263 : }
1264 :
1265 : //No point in reading data if its not matching in frequency
1266 0 : if(max(chanMap)==-1)
1267 0 : return;
1268 0 : Cube<Complex> data;
1269 0 : Cube<Int> flags;
1270 0 : getInterpolateArrays(vb, data, flags);
1271 :
1272 : Complex *datStorage;
1273 : Bool isCopy;
1274 0 : datStorage=data.getStorage(isCopy);
1275 : // NOTE: with MS V2.0 the pointing could change per antenna and timeslot
1276 : //
1277 0 : Vector<Int> rowFlags(vb.flagRow().nelements());
1278 0 : rowFlags=0;
1279 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1280 0 : if(vb.flagRow()(rownr)) rowFlags(rownr)=1;
1281 : //single dish yes ?
1282 0 : 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 0 : Matrix<Double> xyPositions(2, endRow-startRow+1);
1334 0 : xyPositions=-1e9;
1335 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1336 0 : if(getXYPos(vb, rownr)) {
1337 0 : xyPositions(0, rownr)=xyPos(0);
1338 0 : xyPositions(1, rownr)=xyPos(1);
1339 : }
1340 : }
1341 :
1342 : Bool del;
1343 : // IPosition s(data.shape());
1344 0 : const IPosition& fs=data.shape();
1345 0 : std::vector<Int> s(fs.begin(), fs.end());
1346 0 : dgridsd(xyPositions.getStorage(del),
1347 : datStorage,
1348 0 : &s[0],
1349 0 : &s[1],
1350 0 : flags.getStorage(del),
1351 0 : rowFlags.getStorage(del),
1352 0 : &s[2],
1353 : &row,
1354 0 : 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 0 : data.putStorage(datStorage, isCopy);
1366 : }
1367 0 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1368 : }
1369 :
1370 : // Helper functions for SDGrid::makeImage
1371 : namespace {
1372 : inline
1373 0 : void setupVisBufferForFTMachineType(FTMachine::Type type, VisBuffer& vb) {
1374 0 : 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 0 : case FTMachine::COVERAGE:
1383 0 : vb.visCube() = Complex(1.0);
1384 0 : break;
1385 0 : default:
1386 0 : break;
1387 : }
1388 0 : }
1389 :
1390 : inline
1391 0 : 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 0 : auto haveCorrectedData = not (vi.msColumns().correctedData().isNull());
1396 0 : auto tunedType =
1397 0 : ((in_type == FTMachine::CORRECTED) && (not haveCorrectedData)) ?
1398 : FTMachine::OBSERVED : in_type;
1399 :
1400 : // Compute output parameters
1401 0 : switch(tunedType) {
1402 0 : case FTMachine::RESIDUAL:
1403 : case FTMachine::MODEL:
1404 : case FTMachine::CORRECTED:
1405 0 : out_dopsf = false;
1406 0 : out_type = tunedType;
1407 0 : break;
1408 0 : case FTMachine::PSF:
1409 : case FTMachine::COVERAGE:
1410 0 : out_dopsf = true;
1411 0 : out_type = tunedType;
1412 0 : break;
1413 0 : default:
1414 0 : out_dopsf = false;
1415 0 : out_type = FTMachine::OBSERVED;
1416 0 : break;
1417 : }
1418 0 : }
1419 :
1420 0 : void abortOnPolFrameChange(const StokesImageUtil::PolRep refPolRep, const String & refMsName, ROVisibilityIterator &vi) {
1421 0 : auto vb = vi.getVisBuffer();
1422 0 : const auto polRep = (vb->polFrame() == MSIter::Linear) ?
1423 0 : StokesImageUtil::LINEAR : StokesImageUtil::CIRCULAR;
1424 0 : const auto polFrameChanged = (polRep != refPolRep);
1425 0 : 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 0 : }
1461 :
1462 : } // SDGrid::makeImage helper functions namespace
1463 :
1464 0 : void SDGrid::nextChunk(ROVisibilityIterator &vi) {
1465 : #if defined(SDGRID_PERFS)
1466 0 : StartStop trigger(cNextChunk);
1467 : #endif
1468 0 : vi.nextChunk();
1469 0 : }
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 0 : 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 0 : 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 0 : vi.origin();
1487 0 : const auto imgPolRep = (vb.polFrame() == MSIter::Linear) ?
1488 0 : StokesImageUtil::LINEAR : StokesImageUtil::CIRCULAR;
1489 0 : StokesImageUtil::changeCStokesRep(theImage, imgPolRep);
1490 0 : const auto firstMsName = vb.msName();
1491 :
1492 0 : initializeToSky(theImage, weight, vb);
1493 :
1494 : // Setup SDGrid Cache Manager
1495 0 : const auto onDuty = cacheIsEnabled;
1496 0 : const auto accessMode = cache.isEmpty() ? Cache::AccessMode::WRITE
1497 0 : : Cache::AccessMode::READ;
1498 0 : CacheManager cacheManager(cache, onDuty, accessMode);
1499 :
1500 : // Loop over the visibilities, putting VisBuffers
1501 0 : for (vi.originChunks(); vi.moreChunks(); nextChunk(vi)) {
1502 0 : abortOnPolFrameChange(imgPolRep, firstMsName, vi);
1503 : FTMachine::Type actualType;
1504 : Bool doPSF;
1505 0 : if (vi.newMS()) { // Note: the first MS is a new MS
1506 0 : getParamsForFTMachineType(vi, inType, doPSF, actualType);
1507 0 : handleNewMs(vi, theImage);
1508 : }
1509 0 : for (vi.origin(); vi.more(); vi++) {
1510 0 : setupVisBufferForFTMachineType(actualType, vb);
1511 0 : constexpr Int allVbRows = -1;
1512 0 : put(vb, allVbRows, doPSF, actualType);
1513 : }
1514 : }
1515 0 : finalizeToSky();
1516 :
1517 : // Normalize by dividing out weights, etc.
1518 0 : auto doNormalize = (inType != FTMachine::COVERAGE);
1519 0 : getImage(weight, doNormalize);
1520 :
1521 : // Warning message
1522 0 : if (allEQ(weight, 0.0f)) {
1523 0 : LogIO logger(LogOrigin(name(),"makeImage"));
1524 : logger << LogIO::SEVERE
1525 : << "No useful data in SDGrid: all weights are zero"
1526 0 : << LogIO::POST;
1527 : }
1528 0 : }
1529 :
1530 : // Finalize : optionally normalize by weight image
1531 0 : ImageInterface<Complex>& SDGrid::getImage(Matrix<Float>& weights,
1532 : Bool normalize)
1533 : {
1534 0 : AlwaysAssert(lattice, AipsError);
1535 0 : AlwaysAssert(image, AipsError);
1536 :
1537 0 : logIO() << LogOrigin("SDGrid", "getImage") << LogIO::NORMAL;
1538 :
1539 : // execute minmax clipping
1540 0 : clipMinMax();
1541 :
1542 0 : weights.resize(sumWeight.shape());
1543 :
1544 0 : 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 0 : if(normalize) {
1560 0 : 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 0 : 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 0 : IPosition pos(4);
1576 0 : for (Int chan=0; chan < nchan; ++chan){
1577 0 : pos[3]=chan;
1578 0 : for( Int pol=0; pol < npol; ++ pol){
1579 0 : pos[2]=pol;
1580 0 : for (Int y=0; y < ny ; ++y){
1581 0 : pos[1]=y;
1582 0 : for (Int x=0; x < nx; ++x){
1583 0 : pos[0]=x;
1584 0 : Float wgt=wGriddedData(pos);
1585 0 : if(wgt > minWeight_p)
1586 0 : griddedData(pos)=griddedData(pos)/wgt;
1587 : else
1588 0 : 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 0 : IPosition gridShape(4, nx, ny, npol, nchan);
1600 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
1601 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1602 0 : IPosition stride(4, 1);
1603 0 : IPosition trc(blc+image->shape()-stride);
1604 0 : LCBox gridBox(blc, trc, gridShape);
1605 0 : SubLattice<Complex> gridSub(*arrayLattice, gridBox);
1606 :
1607 : // Do the copy
1608 0 : image->copyData(gridSub);
1609 : }
1610 0 : 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 0 : void SDGrid::ok() {
1628 0 : AlwaysAssert(image, AipsError);
1629 0 : }
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 0 : Int SDGrid::getIndex(const MSPointingColumns& mspc, const Double& time,
1640 : const Double& interval, const Int& antid) {
1641 : //Int start=lastIndex_p;
1642 0 : Int start=lastIndexPerAnt_p[antid];
1643 0 : Double tol=interval < 0.0 ? time*C::dbl_epsilon : interval/2.0;
1644 : // Search forwards
1645 0 : Int nrows=mspc.time().nrow();
1646 0 : 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 0 : if (antid >= 0 && mspc.antennaId()(i) != antid) {
1650 0 : continue;
1651 : }
1652 :
1653 0 : 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 0 : Double mspc_interval = mspc.interval()(i);
1658 :
1659 0 : 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 0 : if(abs(midpoint-time) <= (mspc_interval/2.0+tol)) {
1668 : //lastIndex_p=i;
1669 0 : lastIndexPerAnt_p[antid]=i;
1670 0 : return i;
1671 : }
1672 : }
1673 : }
1674 : // Search backwards
1675 0 : for (Int i=start;i>=0;i--) {
1676 0 : if (antid >= 0 && mspc.antennaId()(i) != antid) {
1677 0 : continue;
1678 : }
1679 0 : Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
1680 0 : Double mspc_interval = mspc.interval()(i);
1681 0 : 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 0 : if(abs(midpoint-time) <= (mspc_interval/2.0+tol)) {
1690 : //lastIndex_p=i;
1691 0 : lastIndexPerAnt_p[antid]=i;
1692 0 : return i;
1693 : }
1694 : }
1695 : }
1696 : // No match!
1697 0 : return -1;
1698 : }
1699 :
1700 0 : Bool SDGrid::getXYPos(const VisBuffer& vb, Int row) {
1701 :
1702 : // Cache control
1703 0 : if (cacheIsEnabled and cache.isReadable()) {
1704 0 : cache.loadRowPixel();
1705 0 : return rowPixel.isValid;
1706 : }
1707 0 : const auto onDuty = cacheIsEnabled and cache.isWriteable();
1708 0 : CacheWriter cacheWriter(cache, onDuty);
1709 :
1710 : // Until we manage to compute a valid one ...
1711 0 : rowPixel.isValid = false;
1712 :
1713 : // Select the POINTING table (columns) we'll work with.
1714 0 : const auto haveConvertedColumn = ramPointingTable.nrow() > 0;
1715 0 : const MSPointingColumns& pointingColumns = haveConvertedColumn ?
1716 0 : *ramPointingColumnsPtr
1717 0 : : vb.msColumns().pointing();
1718 :
1719 0 : const auto nPointings = pointingColumns.nrow();
1720 0 : const auto havePointings = (nPointings >= 1);
1721 :
1722 : // We'll need to call these many times, so let's call them once for good
1723 0 : const auto rowTime = vb.time()(row);
1724 0 : const auto rowTimeInterval = vb.timeInterval()(row);
1725 0 : 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 0 : constexpr Int invalidIndex = -1;
1732 0 : Int pointingIndex = invalidIndex;
1733 0 : if (havePointings) {
1734 : #if defined(SDGRID_PERFS)
1735 0 : 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 0 : if (vb.msId() != msId_p) {
1748 0 : lastIndex_p = 0; // no longer used
1749 0 : if (lastIndexPerAnt_p.nelements() < (size_t)vb.numberAnt()) {
1750 0 : lastIndexPerAnt_p.resize(vb.numberAnt());
1751 : }
1752 0 : lastIndexPerAnt_p = 0;
1753 0 : msId_p = vb.msId();
1754 0 : 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 0 : constexpr Double useTinyTolerance = -1.0;
1762 0 : pointingIndex = getIndex(pointingColumns, rowTime, useTinyTolerance , rowAntenna1);
1763 :
1764 0 : auto foundPointing = (pointingIndex >= 0);
1765 0 : 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 0 : const auto foundValidPointing = (foundPointing and (static_cast<rownr_t>(pointingIndex) < nPointings));
1774 :
1775 0 : 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 0 : Bool needInterpolation = False;
1793 0 : if (havePointings) {
1794 0 : const auto pointingInterval = pointingColumns.interval()(pointingIndex);
1795 0 : if (rowTimeInterval < pointingInterval) needInterpolation = True;
1796 : }
1797 0 : const auto mustInterpolate = havePointings && needInterpolation;
1798 :
1799 : // 3. Create interpolator if needed
1800 0 : if (mustInterpolate) {
1801 0 : if (not isSplineInterpolationReady) {
1802 : #if defined(SDGRID_PERFS)
1803 0 : StartStop trigger(cComputeSplines);
1804 : #endif
1805 : const auto nAntennas = static_cast<size_t>(
1806 0 : vb.msColumns().antenna().nrow()
1807 : );
1808 0 : interpolator = new SDPosInterpolator(
1809 : pointingColumns,
1810 0 : pointingDirCol_p,
1811 : nAntennas
1812 0 : );
1813 0 : isSplineInterpolationReady = true;
1814 : } else {
1815 : // We have an interpolator. Re-use it if possible.
1816 0 : const auto canReuseInterpolator = interpolator->inTimeRange(rowTime, rowAntenna1);
1817 0 : if (not canReuseInterpolator) {
1818 : // setup spline interpolator for the current dataset (CAS-11261, 2018/5/22 WK)
1819 : // delete and re-create it
1820 0 : delete interpolator;
1821 0 : interpolator = 0;
1822 : const auto nAntennas = static_cast<size_t>(
1823 0 : vb.msColumns().antenna().nrow()
1824 : );
1825 0 : interpolator = new SDPosInterpolator(
1826 : pointingColumns,
1827 0 : pointingDirCol_p,
1828 : nAntennas
1829 0 : );
1830 : }
1831 : }
1832 : }
1833 :
1834 : // 4. Create the direction conversion machine if needed
1835 :
1836 0 : if ( pointingDirCol_p == "SOURCE_OFFSET" or
1837 0 : 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 0 : const auto needDirectionConverter = (
1844 0 : not havePointings or not haveConvertedColumn or fixMovingSource_p
1845 : );
1846 :
1847 0 : if (needDirectionConverter) {
1848 0 : if (not pointingToImage) {
1849 : // Set the frame
1850 : const auto & rowAntenna1Position =
1851 0 : vb.msColumns().antenna().positionMeas()(rowAntenna1);
1852 : // set dummy time stamp 1 day before rowTime
1853 0 : const MEpoch dummyEpoch(Quantity(rowTime - 86400.0, "s"));
1854 0 : mFrame_p = MeasFrame(dummyEpoch, rowAntenna1Position);
1855 :
1856 : // Remember antenna id for next call,
1857 : // which may be done using a different VisBuffer ...
1858 0 : lastAntID_p = rowAntenna1;
1859 :
1860 : // Compute the "model" required to setup the conversion machine
1861 0 : if (havePointings) {
1862 0 : worldPosMeas = mustInterpolate ?
1863 : directionMeas(pointingColumns, pointingIndex, rowTime)
1864 0 : : 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 0 : MDirection::Ref outRef(directionCoord.directionType(), mFrame_p);
1873 0 : pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
1874 0 : 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 0 : MDirection _dir_tmp = (*pointingToImage)();
1880 : }
1881 : }
1882 :
1883 : // 5. Update the frame holding the measurements for this row
1884 0 : const MEpoch rowEpoch(Quantity(rowTime, "s"));
1885 : // ---- Always reset the epoch
1886 : {
1887 : #if defined(SDGRID_PERFS)
1888 0 : StartStop trigger(cResetFrame);
1889 : #endif
1890 0 : mFrame_p.resetEpoch(rowEpoch);
1891 : }
1892 : // ---- Reset antenna position only if antenna changed since we were last called
1893 0 : if (lastAntID_p != rowAntenna1) {
1894 : // Debug messages
1895 0 : if (lastAntID_p == -1) {
1896 : // antenna ID is unset
1897 0 : logIO_p << LogIO::DEBUGGING
1898 : << "updating antenna position for conversion: new MS ID " << msId_p
1899 0 : << ", 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 0 : vb.msColumns().antenna().positionMeas()(rowAntenna1);
1908 : {
1909 : #if defined(SDGRID_PERFS)
1910 0 : StartStop trigger(cResetFrame);
1911 : #endif
1912 0 : mFrame_p.resetPosition(rowAntenna1Position);
1913 : }
1914 : // Remember antenna id for next call,
1915 : // which may be done using a different VisBuffer ...
1916 0 : lastAntID_p = rowAntenna1;
1917 : }
1918 :
1919 : // 6. Compute user-specified column direction at data-taking time,
1920 : // in image's direction reference frame
1921 0 : if (havePointings) {
1922 0 : if (mustInterpolate) {
1923 : #if defined(SDGRID_PERFS)
1924 0 : cInterpolateDirection.start();
1925 : #endif
1926 : const auto interpolatedDirection =
1927 0 : directionMeas(pointingColumns, pointingIndex, rowTime);
1928 : #if defined(SDGRID_PERFS)
1929 0 : cInterpolateDirection.stop();
1930 0 : cConvertDirection.start();
1931 : #endif
1932 : worldPosMeas = haveConvertedColumn ?
1933 : interpolatedDirection
1934 0 : : (*pointingToImage)(interpolatedDirection);
1935 : #if defined(SDGRID_PERFS)
1936 0 : 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 0 : const auto columnDirection = directionMeas(pointingColumns, pointingIndex);
1947 : worldPosMeas = haveConvertedColumn ?
1948 : columnDirection
1949 0 : : (*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 0 : StartStop trigger(cComputeDirectionPixel);
1960 : #endif
1961 :
1962 0 : rowPixel.isValid = directionCoord.toPixel(xyPos, worldPosMeas);
1963 : }
1964 :
1965 0 : 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 0 : if (fixMovingSource_p) {
1975 : #if defined(SDGRID_PERFS)
1976 0 : StartStop trigger(cHandleMovingSource);
1977 : #endif
1978 0 : if (xyPosMovingOrig_p.nelements() < 2) {
1979 0 : directionCoord.toPixel(xyPosMovingOrig_p, firstMovingDir_p);
1980 : }
1981 : //via HADEC or AZEL for parallax of near sources
1982 0 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
1983 0 : MDirection tmphadec = MDirection::Convert(movingDir_p, outref1)();
1984 0 : MDirection actSourceDir = (*pointingToImage)(tmphadec);
1985 0 : Vector<Double> actPix;
1986 0 : directionCoord.toPixel(actPix, actSourceDir);
1987 :
1988 : //cout << row << " scan " << vb.scan()(row) << "xyPos " << xyPos
1989 : // << " xyposmovorig " << xyPosMovingOrig_p << " actPix " << actPix << endl;
1990 :
1991 0 : xyPos = xyPos + xyPosMovingOrig_p - actPix;
1992 : }
1993 :
1994 0 : return rowPixel.isValid;
1995 : }
1996 :
1997 0 : MDirection SDGrid::directionMeas(const MSPointingColumns& mspc, const Int& index){
1998 0 : if (pointingDirCol_p == "TARGET") {
1999 0 : return mspc.targetMeas(index);
2000 0 : } 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 0 : } 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 0 : } 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 0 : 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 0 : MDirection SDGrid::directionMeas(const MSPointingColumns& mspc, const Int& index,
2025 : const Double& time){
2026 : //spline interpolation
2027 0 : if (isSplineInterpolationReady) {
2028 0 : Int antid = mspc.antennaId()(index);
2029 0 : if (interpolator->doSplineInterpolation(antid)) {
2030 0 : 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 0 : void SDGrid::pickWeights(const VisBuffer& vb, Matrix<Float>& weight){
2136 : #if defined(SDGRID_PERFS)
2137 0 : StartStop trigger(cPickWeights);
2138 : #endif
2139 : //break reference
2140 0 : weight.resize();
2141 :
2142 0 : if (useImagingWeight_p) {
2143 0 : weight.reference(vb.imagingWeight());
2144 : } else {
2145 0 : const Cube<Float> weightspec(vb.weightSpectrum());
2146 0 : weight.resize(vb.nChannel(), vb.nRow());
2147 :
2148 0 : if (weightspec.nelements() == 0) {
2149 0 : for (Int k = 0; k < vb.nRow(); ++k) {
2150 : //cerr << "nrow " << vb.nRow() << " " << weight.shape() << " " << weight.column(k).shape() << endl;
2151 0 : 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 0 : }
2171 :
2172 0 : void SDGrid::clipMinMax() {
2173 0 : if (clipminmax_) {
2174 : Bool gmin_b, gmax_b, wmin_b, wmax_b, np_b;
2175 0 : const auto *gmin_p = gmin_.getStorage(gmin_b);
2176 0 : const auto *gmax_p = gmax_.getStorage(gmax_b);
2177 0 : const auto *wmin_p = wmin_.getStorage(wmin_b);
2178 0 : const auto *wmax_p = wmax_.getStorage(wmax_b);
2179 0 : const auto *np_p = npoints_.getStorage(np_b);
2180 :
2181 : Bool data_b, weight_b, sumw_b;
2182 0 : auto data_p = griddedData.getStorage(data_b);
2183 0 : auto weight_p = wGriddedData.getStorage(weight_b);
2184 0 : auto sumw_p = sumWeight.getStorage(sumw_b);
2185 :
2186 0 : auto arrayShape = griddedData.shape();
2187 0 : size_t num_xy = arrayShape.getFirst(2).product();
2188 0 : size_t num_polchan = arrayShape.getLast(2).product();
2189 0 : for (size_t i = 0; i < num_xy; ++i) {
2190 0 : for (size_t j = 0; j < num_polchan; ++j) {
2191 0 : auto k = i * num_polchan + j;
2192 0 : if (np_p[k] == 1) {
2193 0 : auto wt = wmin_p[k];
2194 0 : data_p[k] = wt * gmin_p[k];
2195 0 : weight_p[k] = wt;
2196 0 : sumw_p[j] += wt;
2197 0 : } else if (np_p[k] == 2) {
2198 0 : auto wt = wmin_p[k];
2199 0 : data_p[k] = wt * gmin_p[k];
2200 0 : weight_p[k] = wt;
2201 0 : sumw_p[j] += wt;
2202 0 : wt = wmax_p[k];
2203 0 : data_p[k] += wt * gmax_p[k];
2204 0 : weight_p[k] += wt;
2205 0 : sumw_p[j] += wt;
2206 : }
2207 : }
2208 : }
2209 :
2210 0 : wGriddedData.putStorage(weight_p, weight_b);
2211 0 : griddedData.putStorage(data_p, data_b);
2212 0 : sumWeight.putStorage(sumw_p, sumw_b);
2213 :
2214 0 : npoints_.freeStorage(np_p, np_b);
2215 0 : wmax_.freeStorage(wmax_p, wmax_b);
2216 0 : wmin_.freeStorage(wmin_p, wmin_b);
2217 0 : gmax_.freeStorage(gmax_p, gmax_b);
2218 0 : gmin_.freeStorage(gmin_p, gmin_b);
2219 : }
2220 0 : }
2221 :
2222 0 : SDGrid::MaskedPixelRef::MaskedPixelRef(Vector<Double>& xyIn, Bool isValidIn)
2223 : : xy {xyIn},
2224 0 : isValid {isValidIn}
2225 0 : {}
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 0 : SDGrid::MaskedPixel::MaskedPixel(Double xIn, Double yIn, Bool isValidIn)
2235 : : x {xIn},
2236 : y {yIn},
2237 0 : isValid {isValidIn}
2238 0 : {}
2239 :
2240 0 : SDGrid::Cache::Cache(SDGrid &parent)
2241 : : sdgrid {parent},
2242 : isOpened {false},
2243 : accessMode {AccessMode::READ},
2244 : canRead {false},
2245 : canWrite {false},
2246 0 : inputPixel {sdgrid.rowPixel},
2247 0 : outputPixel {sdgrid.rowPixel}
2248 0 : {}
2249 :
2250 0 : const casacore::String& SDGrid::Cache::className() {
2251 0 : static casacore::String className_ = "SDGrid::Cache";
2252 0 : 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 0 : SDGrid::Cache::open(AccessMode accessModeIn) {
2274 0 : 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 0 : isOpened = true;
2279 0 : accessMode = accessModeIn;
2280 0 : canRead = accessMode == AccessMode::READ;
2281 0 : canWrite = accessMode == AccessMode::WRITE;
2282 0 : if (outputPixel.xy.size() < 2) {
2283 0 : outputPixel.xy.resize(2);
2284 0 : outputPixel.isValid = false;
2285 : }
2286 0 : rewind();
2287 0 : }
2288 :
2289 : void
2290 0 : SDGrid::Cache::rewind() {
2291 : // Writer
2292 0 : msPixels = nullptr;
2293 :
2294 : // Reader
2295 0 : msCacheReadIterator = msCaches.cbegin();
2296 0 : pixelReadIterator = Pixels::const_iterator();
2297 0 : }
2298 :
2299 : void
2300 0 : SDGrid::Cache::close() {
2301 0 : 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 0 : if (isWriteable()) {
2306 : // Make sure we have 1 pixel per row
2307 0 : for (const auto & msCache : msCaches) {
2308 0 : 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 0 : isOpened = false;
2323 0 : accessMode = AccessMode::READ;
2324 0 : canRead = false;
2325 0 : canWrite = false;
2326 0 : }
2327 :
2328 : Bool
2329 0 : SDGrid::Cache::isReadable() const {
2330 0 : return isOpened and canRead;
2331 : }
2332 :
2333 : Bool
2334 0 : SDGrid::Cache::isWriteable() const {
2335 0 : return isOpened and canWrite;
2336 : }
2337 :
2338 : Bool
2339 0 : SDGrid::Cache::isEmpty() const {
2340 0 : return msCaches.empty();
2341 : }
2342 :
2343 : void
2344 0 : SDGrid::Cache::clear() {
2345 0 : msCaches.clear();
2346 0 : }
2347 :
2348 0 : SDGrid::Cache::MsCache::MsCache(const String& msPathIn, const String& msTableNameIn, rownr_t nRowsIn)
2349 : : msPath {msPathIn},
2350 : msTableName {msTableNameIn},
2351 0 : nRows {nRowsIn}
2352 : {
2353 0 : pixels.reserve(nRows);
2354 0 : }
2355 :
2356 0 : Bool SDGrid::Cache::MsCache::isConsistent() const {
2357 0 : return pixels.size() == nRows;
2358 : }
2359 :
2360 : void
2361 0 : SDGrid::Cache::newMS(const MeasurementSet& ms) {
2362 0 : LogIO os(LogOrigin(className(),"newMS"));
2363 0 : const auto msPath = ms.getPartNames()[0];
2364 :
2365 0 : if (isWriteable()) {
2366 0 : os << "Will compute and cache spectra pixels coordinates for: " << msPath << LogIO::POST;
2367 0 : msCaches.emplace_back(msPath, ms.tableName(), ms.nrow());
2368 0 : msPixels = &(msCaches.back().pixels);
2369 0 : return;
2370 : }
2371 :
2372 0 : if (isReadable()) {
2373 0 : os << "Will load cached spectra pixels coordinates for: " << msPath << LogIO::POST;
2374 0 : if (msCacheReadIterator == msCaches.cend()) {
2375 0 : os << "BUG! Cached data missing for: " << msPath << LogIO::EXCEPTION;
2376 : }
2377 0 : const auto & pixelsToLoad = msCacheReadIterator->pixels;
2378 0 : 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 0 : 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 0 : pixelReadIterator = pixelsToLoad.cbegin();
2389 0 : ++msCacheReadIterator;
2390 : }
2391 : }
2392 :
2393 : void
2394 0 : SDGrid::Cache::storeRowPixel() {
2395 0 : msPixels->emplace_back(
2396 0 : inputPixel.xy[0],
2397 0 : inputPixel.xy[1],
2398 0 : inputPixel.isValid
2399 : );
2400 0 : }
2401 :
2402 : void
2403 0 : SDGrid::Cache::loadRowPixel() {
2404 0 : const auto & pixel = *pixelReadIterator;
2405 0 : outputPixel.xy[0] = pixel.x;
2406 0 : outputPixel.xy[1] = pixel.y;
2407 0 : outputPixel.isValid = pixel.isValid;
2408 0 : ++pixelReadIterator;
2409 0 : }
2410 :
2411 0 : SDGrid::CacheManager::CacheManager(Cache& cacheIn,
2412 0 : Bool onDutyIn, Cache::AccessMode accessModeIn)
2413 : : cache {cacheIn},
2414 : onDuty {onDutyIn},
2415 0 : accessMode {accessModeIn}
2416 : {
2417 0 : if (onDuty) {
2418 0 : cache.open(accessMode);
2419 : }
2420 0 : }
2421 :
2422 0 : SDGrid::CacheManager::~CacheManager()
2423 : {
2424 0 : if (onDuty) {
2425 0 : cache.close();
2426 : }
2427 0 : }
2428 :
2429 0 : SDGrid::CacheWriter::CacheWriter(Cache& cacheIn,
2430 0 : Bool onDutyIn)
2431 : : cache {cacheIn},
2432 0 : onDuty {onDutyIn}
2433 0 : {}
2434 :
2435 0 : SDGrid::CacheWriter::~CacheWriter()
2436 : {
2437 0 : if (onDuty) {
2438 0 : cache.storeRowPixel();
2439 : }
2440 0 : }
2441 :
2442 0 : const String & SDGrid::toString(const ConvertFirst convertFirst) {
2443 : static const std::array<String,3> name {
2444 : "NEVER",
2445 : "ALWAYS",
2446 : "AUTO"
2447 0 : };
2448 0 : switch(convertFirst){
2449 0 : case ConvertFirst::NEVER:
2450 : case ConvertFirst::ALWAYS:
2451 : case ConvertFirst::AUTO:
2452 0 : 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 0 : 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 0 : for ( const auto scheme : schemes ) {
2474 0 : 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 0 : void SDGrid::setConvertFirst(const casacore::String &name) {
2489 0 : processingScheme = fromString(name);
2490 0 : }
2491 :
2492 0 : Bool SDGrid::mustConvertPointingColumn(
2493 : const MeasurementSet &ms
2494 : ) {
2495 : const auto haveCachedSpectraPixelCoordinates =
2496 0 : cacheIsEnabled and cache.isReadable();
2497 0 : if (haveCachedSpectraPixelCoordinates) return False;
2498 :
2499 0 : switch(processingScheme){
2500 0 : case ConvertFirst::ALWAYS: return True;
2501 0 : 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 0 : void SDGrid::handleNewMs(
2525 : ROVisibilityIterator &vi,
2526 : const ImageInterface<Complex>& image) {
2527 :
2528 : // Synchronize spatial coordinates cache
2529 0 : if (cacheIsEnabled) cache.newMS(vi.ms());
2530 :
2531 : // Handle interpolate-convert processing scheme
2532 0 : 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 0 : ramPointingTable = MSPointing();
2540 0 : ramPointingColumnsPtr.reset();
2541 : }
2542 0 : }
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
|
