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