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 0 : ChronoStat::ChronoStat(const string & name)
93 0 : : name_ {name},
94 0 : started_ {false},
95 0 : n_laps_ {0},
96 0 : n_overflows_ {0},
97 0 : n_underflows_ {0},
98 0 : laps_sum_ {Duration::zero()},
99 0 : laps_min_ {Duration::max()},
100 0 : laps_max_ {Duration::min()}
101 0 : {}
102 :
103 0 : const std::string& ChronoStat::name() const {
104 0 : return name_;
105 : }
106 :
107 0 : void ChronoStat::setName(const std::string& name) {
108 0 : name_ = name;
109 0 : }
110 :
111 0 : void ChronoStat::start() {
112 0 : if (not started_) {
113 0 : started_ = true;
114 0 : ++n_laps_;
115 0 : lap_start_time_ = Clock::now();
116 : } else {
117 0 : ++n_overflows_;
118 : }
119 0 : }
120 0 : void ChronoStat::stop() {
121 0 : if (started_) {
122 0 : auto lap_duration = Clock::now() - lap_start_time_;
123 0 : laps_sum_ += lap_duration;
124 0 : started_ = false;
125 0 : if (lap_duration < laps_min_) laps_min_ = lap_duration;
126 0 : if (lap_duration > laps_max_) laps_max_ = lap_duration;
127 : } else {
128 0 : ++n_underflows_;
129 : }
130 0 : }
131 :
132 0 : ChronoStat::Duration ChronoStat::lapsSum() const {
133 0 : return laps_sum_;
134 : }
135 :
136 0 : ChronoStat::Duration ChronoStat::lapsMin() const {
137 0 : return n_laps_ > 0 ? laps_min_ : Duration::zero();
138 : }
139 :
140 0 : ChronoStat::Duration ChronoStat::lapsMax() const {
141 0 : return n_laps_ > 0 ? laps_max_ : Duration::zero();
142 : }
143 :
144 0 : ChronoStat::Duration ChronoStat::lapsMean() const {
145 0 : return n_laps_ > 0 ? laps_sum_ / n_laps_ : Duration::zero();
146 : }
147 :
148 0 : unsigned int ChronoStat::lapsCount() const {
149 0 : return n_laps_;
150 : }
151 :
152 0 : bool ChronoStat::isEmpty() const {
153 0 : return n_laps_ == 0;
154 : }
155 :
156 0 : unsigned int ChronoStat::nOverflows() const {
157 0 : return n_overflows_;
158 : }
159 :
160 0 : unsigned int ChronoStat::nUnderflows() const {
161 0 : return n_underflows_;
162 : }
163 :
164 0 : std::string ChronoStat::quote(const std::string& s) const {
165 0 : return "\"" + s + "\"";
166 : }
167 :
168 0 : std::string ChronoStat::json() const {
169 0 : std::ostringstream os;
170 0 : os << quote(name()) << ": {"
171 0 : << quote("sum") << ": " << lapsSum().count()
172 0 : << " ," << quote("count") << ": " << lapsCount()
173 0 : << " ," << quote("min") << ": " << lapsMin().count()
174 0 : << " ," << quote("mean") << ": " << lapsMean().count()
175 0 : << " ," << quote("max") << ": " << lapsMax().count()
176 0 : << " ," << quote("overflows") << ": " << nOverflows()
177 0 : << " ," << quote("underflows") << ": " << nUnderflows()
178 0 : << "}";
179 0 : return os.str();
180 0 : }
181 :
182 0 : std::ostream& operator<<(std::ostream &os, const ChronoStat &c) {
183 0 : constexpr auto eol = '\n';
184 0 : os << "name: " << c.name() << eol
185 0 : << "sum: " << std::setw(20) << std::right << c.lapsSum().count() << eol
186 0 : << "count: " << std::setw(20) << std::right << c.lapsCount() << eol
187 0 : << "min: " << std::setw(20) << std::right << c.lapsMin().count() << eol
188 0 : << "mean: " << std::setw(20) << std::right << c.lapsMean().count() << eol
189 0 : << "max: " << std::setw(20) << std::right << c.lapsMax().count() << eol
190 0 : << "overflows: " << std::setw(20) << std::right << c.nOverflows() << eol
191 0 : << "underflows: " << std::setw(20) << std::right << c.nUnderflows();
192 0 : return os;
193 : }
194 :
195 0 : StartStop::StartStop(ChronoStat &c)
196 0 : : c_ {c} {
197 0 : c_.start();
198 0 : }
199 :
200 0 : StartStop::~StartStop() {
201 0 : c_.stop();
202 0 : }
203 :
204 : } // namespace sdgrid_perfs
205 :
206 : using namespace sdgrid_perfs;
207 :
208 : #endif
209 :
210 0 : SDGrid::SDGrid(SkyJones& sj, Int icachesize, Int itilesize,
211 0 : String iconvType, Int userSupport, Bool useImagingWeight)
212 0 : : FTMachine(), sj_p(&sj), imageCache(0), wImageCache(0),
213 0 : cachesize(icachesize), tilesize(itilesize),
214 0 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
215 0 : pointingToImage(0), rowPixel {xyPos, false}, userSetSupport_p(userSupport),
216 0 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
217 0 : minWeight_p(0.), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
218 0 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(false),
219 0 : cache {Cache(*(const_cast<SDGrid *>(this)))},
220 0 : cacheIsEnabled {false}
221 : {
222 0 : lastIndex_p=0;
223 0 : initPerfs();
224 0 : }
225 :
226 0 : SDGrid::SDGrid(MPosition& mLocation, SkyJones& sj, Int icachesize, Int itilesize,
227 0 : String iconvType, Int userSupport, Float minweight, Bool clipminmax, Bool useImagingWeight)
228 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 0 : void SDGrid::initPerfs() {
298 : #if defined(SDGRID_PERFS)
299 0 : cNextChunk.setName("iterateNextChunk");
300 0 : cMatchAllSpwChans.setName("matchAllSpwChans");
301 0 : cMatchChannel.setName("matchChannel");
302 0 : cPickWeights.setName("pickWeights");
303 0 : cInterpolateFrequencyToGrid.setName("interpolateFrequencyToGrid");
304 0 : cSearchValidPointing.setName("searchValidPointing");
305 0 : cComputeSplines.setName("computeSplines");
306 0 : cResetFrame.setName("resetFrame");
307 0 : cInterpolateDirection.setName("interpolateDirection");
308 0 : cConvertDirection.setName("convertDirection");
309 0 : cComputeDirectionPixel.setName("computeDirectionPixel");
310 0 : cHandleMovingSource.setName("handleMovingSource");
311 0 : cGridData.setName("gridData");
312 : #endif
313 0 : }
314 :
315 :
316 :
317 : //----------------------------------------------------------------------
318 0 : SDGrid& SDGrid::operator=(const SDGrid& other)
319 : {
320 0 : if(this!=&other) {
321 : //Do the base parameters
322 0 : FTMachine::operator=(other);
323 0 : sj_p=other.sj_p;
324 0 : imageCache=other.imageCache;
325 0 : wImage=other.wImage;
326 0 : wImageCache=other.wImageCache;
327 0 : cachesize=other.cachesize;
328 0 : tilesize=other.tilesize;
329 0 : isTiled=other.isTiled;
330 0 : lattice=other.lattice;
331 0 : arrayLattice=other.arrayLattice;
332 0 : wLattice=other.wLattice;
333 0 : wArrayLattice=other.wArrayLattice;
334 0 : convType=other.convType;
335 0 : if(other.wImage !=0)
336 0 : wImage=(TempImage<Float> *)other.wImage->cloneII();
337 : else
338 0 : wImage=0;
339 0 : pointingDirCol_p=other.pointingDirCol_p;
340 0 : pointingToImage=0;
341 0 : xyPos.resize();
342 0 : xyPos=other.xyPos;
343 0 : rowPixel = MaskedPixelRef(xyPos, false);
344 0 : xyPosMovingOrig_p=other.xyPosMovingOrig_p;
345 0 : convFunc.resize();
346 0 : convFunc=other.convFunc;
347 0 : convSampling=other.convSampling;
348 0 : convSize=other.convSize;
349 0 : convSupport=other.convSupport;
350 0 : userSetSupport_p=other.userSetSupport_p;
351 0 : lastIndex_p=0;
352 0 : lastIndexPerAnt_p=0;
353 0 : lastAntID_p=-1;
354 0 : msId_p=-1;
355 0 : useImagingWeight_p=other.useImagingWeight_p;
356 0 : clipminmax_=other.clipminmax_;
357 0 : cacheIsEnabled=false;
358 0 : cache=Cache(*(const_cast<SDGrid *>(this)));
359 : };
360 0 : return *this;
361 : };
362 :
363 0 : String SDGrid::name() const{
364 0 : return String("SDGrid");
365 : }
366 :
367 : void
368 0 : SDGrid::setEnableCache(Bool doEnable) {
369 0 : cacheIsEnabled = doEnable;
370 0 : }
371 :
372 : //----------------------------------------------------------------------
373 : // Odds are that it changed.....
374 0 : Bool SDGrid::changed(const VisBuffer& /*vb*/) {
375 0 : return false;
376 : }
377 :
378 : //----------------------------------------------------------------------
379 0 : SDGrid::SDGrid(const SDGrid& other)
380 : : FTMachine(),
381 0 : rowPixel {MaskedPixelRef(xyPos)},
382 0 : cache {Cache(*const_cast<SDGrid *>(this))}
383 : {
384 0 : operator=(other);
385 0 : }
386 :
387 : #define NEED_UNDERSCORES
388 : #if defined(NEED_UNDERSCORES)
389 : #define grdsf grdsf_
390 : #define grdgauss grdgauss_
391 : #define grdjinc1 grdjinc1_
392 : #endif
393 :
394 : extern "C" {
395 : void grdsf(Double*, Double*);
396 : void grdgauss(Double*, Double*, Double*);
397 : void grdjinc1(Double*, Double*, Int*, Double*);
398 : }
399 :
400 0 : void SDGrid::dumpConvolutionFunction(const casacore::String &outfile,
401 : const casacore::Vector<casacore::Float> &f) const {
402 :
403 0 : ofstream ofs(outfile.c_str());
404 0 : for (size_t i = 0 ; i < f.nelements() ; i++) {
405 0 : ofs << i << " " << f[i] << endl;
406 : }
407 0 : ofs.close();
408 0 : }
409 :
410 : //----------------------------------------------------------------------
411 0 : void SDGrid::init() {
412 :
413 0 : logIO() << LogOrigin("SDGrid", "init") << LogIO::NORMAL;
414 :
415 0 : ok();
416 :
417 0 : isTiled = false;
418 0 : nx = image->shape()(0);
419 0 : ny = image->shape()(1);
420 0 : npol = image->shape()(2);
421 0 : nchan = image->shape()(3);
422 :
423 0 : sumWeight.resize(npol, nchan);
424 :
425 : // Set up image cache needed for gridding. For BOX-car convolution
426 : // we can use non-overlapped tiles. Otherwise we need to use
427 : // overlapped tiles and additive gridding so that only increments
428 : // to a tile are written.
429 0 : if (imageCache) delete imageCache;
430 0 : imageCache = nullptr;
431 :
432 0 : convType = downcase(convType);
433 0 : logIO() << "Convolution function : " << convType << LogIO::DEBUG1 << LogIO::POST;
434 :
435 : // Compute convolution function
436 0 : if (convType=="pb") {
437 : // Do nothing
438 : }
439 0 : else if (convType=="box") {
440 0 : convSupport=(userSetSupport_p >= 0) ? userSetSupport_p : 0;
441 0 : logIO() << "Support : " << convSupport << " pixels" << LogIO::POST;
442 0 : convSampling=100;
443 0 : convSize=convSampling*(2*convSupport+2);
444 0 : convFunc.resize(convSize);
445 0 : convFunc=0.0;
446 0 : for (Int i=0;i<convSize/2;i++) {
447 0 : convFunc(i)=1.0;
448 : }
449 : }
450 0 : else if (convType=="sf") {
451 0 : convSupport=(userSetSupport_p >= 0) ? userSetSupport_p : 3;
452 0 : logIO() << "Support : " << convSupport << " pixels" << LogIO::POST;
453 0 : convSampling=100;
454 0 : convSize=convSampling*(2*convSupport+2);
455 0 : convFunc.resize(convSize);
456 0 : convFunc=0.0;
457 0 : for (Int i=0;i<convSampling*convSupport;i++) {
458 0 : Double nu=Double(i)/Double(convSupport*convSampling);
459 : Double val;
460 0 : grdsf(&nu, &val);
461 0 : convFunc(i)=(1.0-nu*nu)*val;
462 : }
463 : }
464 0 : else if (convType=="gauss") {
465 : // default is b=1.0 (Mangum et al. 2007)
466 0 : Double hwhm=(gwidth_p > 0.0) ? Double(gwidth_p) : sqrt(log(2.0));
467 0 : Float truncate=(truncate_p >= 0.0) ? truncate_p : 3.0*hwhm;
468 0 : convSampling=100;
469 0 : Int itruncate=(Int)(truncate*Double(convSampling)+0.5);
470 0 : logIO() << LogIO::DEBUG1 << "hwhm=" << hwhm << LogIO::POST;
471 0 : logIO() << LogIO::DEBUG1 << "itruncate=" << itruncate << LogIO::POST;
472 : //convSupport=(Int)(truncate+0.5);
473 0 : convSupport = (Int)(truncate);
474 0 : convSupport += (((truncate-(Float)convSupport) > 0.0) ? 1 : 0);
475 0 : convSize=convSampling*(2*convSupport+2);
476 0 : convFunc.resize(convSize);
477 0 : convFunc=0.0;
478 : Double val, x;
479 0 : for (Int i = 0 ; i <= itruncate ; i++) {
480 0 : x = Double(i)/Double(convSampling);
481 0 : grdgauss(&hwhm, &x, &val);
482 0 : convFunc(i)=val;
483 : }
484 : }
485 0 : else if (convType=="gjinc") {
486 : // default is b=2.52, c=1.55 (Mangum et al. 2007)
487 0 : Double hwhm = (gwidth_p > 0.0) ? Double(gwidth_p) : sqrt(log(2.0))*2.52;
488 0 : Double c = (jwidth_p > 0.0) ? Double(jwidth_p) : 1.55;
489 : //Float truncate = truncate_p;
490 0 : convSampling = 100;
491 0 : Int itruncate=(Int)(truncate_p*Double(convSampling)+0.5);
492 0 : logIO() << LogIO::DEBUG1 << "hwhm=" << hwhm << LogIO::POST;
493 0 : logIO() << LogIO::DEBUG1 << "c=" << c << LogIO::POST;
494 0 : logIO() << LogIO::DEBUG1 << "itruncate=" << itruncate << LogIO::POST;
495 : //convSupport=(truncate_p >= 0.0) ? (Int)(truncate_p+0.5) : (Int)(2*c+0.5);
496 0 : Float convSupportF = (truncate_p >= 0.0) ? truncate_p : (2*c);
497 0 : convSupport = (Int)convSupportF;
498 0 : convSupport += (((convSupportF-(Float)convSupport) > 0.0) ? 1 : 0);
499 0 : convSize=convSampling*(2*convSupport+2);
500 0 : convFunc.resize(convSize);
501 0 : convFunc=0.0;
502 : //UNUSED: Double r;
503 : Double x, val1, val2;
504 0 : Int normalize = 1;
505 0 : if (itruncate >= 0) {
506 0 : for (Int i = 0 ; i < itruncate ; i++) {
507 0 : x = Double(i) / Double(convSampling);
508 : //r = Double(i) / (Double(hwhm)*Double(convSampling));
509 0 : grdgauss(&hwhm, &x, &val1);
510 0 : grdjinc1(&c, &x, &normalize, &val2);
511 0 : convFunc(i) = val1 * val2;
512 : }
513 : }
514 : else {
515 : // default is to truncate at first null
516 0 : for (Int i=0;i<convSize;i++) {
517 0 : x = Double(i) / Double(convSampling);
518 : //r = Double(i) / (Double(hwhm)*Double(convSampling));
519 0 : grdjinc1(&c, &x, &normalize, &val2);
520 0 : if (val2 <= 0.0) {
521 : logIO() << LogIO::DEBUG1
522 : << "convFunc is automatically truncated at radius "
523 0 : << x << LogIO::POST;
524 0 : break;
525 : }
526 0 : grdgauss(&hwhm, &x, &val1);
527 0 : convFunc(i) = val1 * val2;
528 : }
529 : }
530 : }
531 : else {
532 0 : logIO_p << "Unknown convolution function " << convType << LogIO::EXCEPTION;
533 : }
534 :
535 : // Convolution function debug
536 : // String outfile = convType + ".dat";
537 : // dumpConvolutionFunction(outfile,convFunc);
538 :
539 0 : if (wImage) delete wImage;
540 0 : wImage = new TempImage<Float>(image->shape(), image->coordinates());
541 :
542 0 : }
543 :
544 0 : void SDGrid::collectPerfs(){
545 : #if defined(SDGRID_PERFS)
546 0 : LogIO os(LogOrigin(name(), "collectPerfs"));
547 : std::vector<std::string> probes_json {
548 : cNextChunk.json()
549 : , cMatchAllSpwChans.json()
550 : , cMatchChannel.json()
551 : , cPickWeights.json()
552 : , cInterpolateFrequencyToGrid.json()
553 : , cSearchValidPointing.json()
554 : , cComputeSplines.json()
555 : , cResetFrame.json()
556 : , cInterpolateDirection.json()
557 : , cConvertDirection.json()
558 : , cComputeDirectionPixel.json()
559 : , cHandleMovingSource.json()
560 : , cGridData.json()
561 0 : };
562 : os << "PERFS<SDGRID> "
563 : << "{ \"note\": \"sum, min, mean, max are in units of nanoseconds.\""
564 : << ", \"probes\": "
565 0 : << "{ " << join(probes_json.data(), probes_json.size(), ", " ) << " }"
566 : << " }"
567 0 : << LogIO::POST;
568 : #endif
569 0 : }
570 :
571 :
572 : // This is nasty, we should use CountedPointers here.
573 0 : SDGrid::~SDGrid() {
574 : //fclose(pfile);
575 0 : if (imageCache) {
576 0 : delete imageCache;
577 0 : imageCache = nullptr;
578 : }
579 0 : if (arrayLattice) {
580 0 : delete arrayLattice;
581 0 : arrayLattice = nullptr;
582 : }
583 0 : if (wImage) {
584 0 : delete wImage;
585 0 : wImage = nullptr;
586 : }
587 0 : if (wImageCache) {
588 0 : delete wImageCache;
589 0 : wImageCache = nullptr;
590 : }
591 0 : if (wArrayLattice) {
592 0 : delete wArrayLattice;
593 0 : wArrayLattice = nullptr;
594 : }
595 0 : if (interpolator) {
596 0 : delete interpolator;
597 0 : interpolator = nullptr;
598 : }
599 :
600 0 : collectPerfs();
601 0 : }
602 :
603 0 : 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 0 : 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 0 : convSupport=max(128, sj_p->support(vb, coords));
613 :
614 0 : convSampling=100;
615 0 : 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 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
621 0 : AlwaysAssert(directionIndex>=0, AipsError);
622 0 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
623 0 : Vector<Double> sampling;
624 0 : sampling = dc.increment();
625 0 : sampling*=1.0/Double(convSampling);
626 0 : 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 0 : uInt row=0;
632 :
633 : // reset lastAntID_p to use correct antenna position
634 0 : lastAntID_p = -1;
635 :
636 0 : const MSPointingColumns& act_mspc = vb.msColumns().pointing();
637 0 : Bool nullPointingTable=(act_mspc.nrow() < 1);
638 : // uInt pointIndex=getIndex(*mspc, vb.time()(row), vb.timeInterval()(row), vb.antenna1()(row));
639 0 : Int pointIndex=-1;
640 0 : if(!nullPointingTable){
641 : //if(vb.newMS()) This thing is buggy...using msId change
642 0 : if(vb.msId() != msId_p){
643 0 : lastIndex_p=0;
644 0 : if(lastIndexPerAnt_p.nelements() < (size_t)vb.numberAnt()) {
645 0 : lastIndexPerAnt_p.resize(vb.numberAnt());
646 : }
647 0 : lastIndexPerAnt_p=0;
648 0 : msId_p=vb.msId();
649 : }
650 0 : pointIndex=getIndex(act_mspc, vb.time()(row), -1.0, vb.antenna1()(row));
651 0 : if(pointIndex < 0)
652 0 : pointIndex=getIndex(act_mspc, vb.time()(row), vb.timeInterval()(row), vb.antenna1()(row));
653 :
654 : }
655 0 : 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 0 : MEpoch epoch(Quantity(vb.time()(row), "s"));
662 0 : if(!pointingToImage) {
663 0 : lastAntID_p=vb.antenna1()(row);
664 0 : MPosition pos;
665 0 : pos=vb.msColumns().antenna().positionMeas()(lastAntID_p);
666 0 : mFrame_p=MeasFrame(epoch, pos);
667 0 : if(!nullPointingTable){
668 0 : 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 0 : MDirection::Ref outRef(dc.directionType(), mFrame_p);
676 0 : pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
677 :
678 0 : if(!pointingToImage) {
679 0 : logIO_p << "Cannot make direction conversion machine" << LogIO::EXCEPTION;
680 : }
681 0 : }
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 0 : if(!nullPointingTable){
696 0 : worldPosMeas=(*pointingToImage)(directionMeas(act_mspc,pointIndex));
697 : }
698 : else{
699 0 : worldPosMeas=(*pointingToImage)(vb.direction1()(row));
700 : }
701 0 : delete pointingToImage;
702 0 : pointingToImage=0;
703 0 : }
704 :
705 0 : directionCoord=coords.directionCoordinate(directionIndex);
706 : //make sure we use the same units
707 0 : 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 0 : Vector<Double> const originalReferencePixel = dc.referencePixel();
715 0 : dc.setReferenceValue(worldPosMeas.getAngle().getValue());
716 : //Vector<Double> unitVec(2);
717 : //unitVec=0.0;
718 : //dc.setReferencePixel(unitVec);
719 0 : Vector<Double> updatedReferencePixel = dc.referencePixel() - originalReferencePixel;
720 0 : dc.setReferencePixel(updatedReferencePixel);
721 :
722 0 : coords.replaceCoordinate(dc, directionIndex);
723 :
724 0 : IPosition pbShape(4, convSize, 2, 1, 1);
725 0 : IPosition start(4, 0, 0, 0, 0);
726 :
727 0 : TempImage<Complex> onedPB(pbShape, coords);
728 :
729 0 : onedPB.set(Complex(1.0, 0.0));
730 :
731 0 : AlwaysAssert(sj_p, AipsError);
732 0 : sj_p->apply(onedPB, onedPB, vb, 0);
733 :
734 0 : IPosition pbSlice(4, convSize, 1, 1, 1);
735 0 : Vector<Float> tempConvFunc=real(onedPB.getSlice(start, pbSlice, true));
736 : // Find number of significant points
737 0 : uInt cfLen=0;
738 0 : for(uInt i=0;i<tempConvFunc.nelements();++i) {
739 0 : if(tempConvFunc(i)<1e-3) break;
740 0 : ++cfLen;
741 : }
742 0 : 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 0 : Vector<Float> trimConvFunc=tempConvFunc(Slice(0,cfLen-1,1));
749 :
750 : // Now fill in the convolution function vector
751 0 : convSupport=cfLen/convSampling;
752 0 : convSize=convSampling*(2*convSupport+2);
753 0 : convFunc.resize(convSize);
754 0 : convFunc=0.0;
755 0 : convFunc(Slice(0,cfLen-1,1))=trimConvFunc(Slice(0,cfLen-1,1));
756 :
757 :
758 0 : }
759 :
760 : // Initialize for a transform from the Sky domain. This means that
761 : // we grid-correct, and FFT the image
762 0 : void SDGrid::initializeToVis(ImageInterface<Complex>& iimage,
763 : const VisBuffer& vb)
764 : {
765 0 : image=&iimage;
766 :
767 0 : ok();
768 :
769 0 : init();
770 :
771 0 : if(convType=="pb") {
772 0 : 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 0 : msId_p = -1;
778 0 : lastAntID_p = -1;
779 :
780 : // Initialize the maps for polarization and channel. These maps
781 : // translate visibility indices into image indices
782 0 : initMaps(vb);
783 :
784 : // First get the CoordinateSystem for the image and then find
785 : // the DirectionCoordinate
786 0 : CoordinateSystem coords=image->coordinates();
787 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
788 0 : AlwaysAssert(directionIndex>=0, AipsError);
789 0 : directionCoord=coords.directionCoordinate(directionIndex);
790 : /*if((image->shape().product())>cachesize) {
791 : isTiled=true;
792 : }
793 : else {
794 : isTiled=false;
795 : }*/
796 0 : isTiled=false;
797 0 : nx = image->shape()(0);
798 0 : ny = image->shape()(1);
799 0 : npol = image->shape()(2);
800 0 : 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 0 : IPosition gridShape(4, nx, ny, npol, nchan);
812 0 : griddedData.resize(gridShape);
813 0 : griddedData = Complex(0.0);
814 :
815 0 : wGriddedData.resize(gridShape);
816 0 : wGriddedData = 0.0;
817 :
818 0 : if(arrayLattice) delete arrayLattice; arrayLattice=0;
819 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
820 :
821 0 : if(wArrayLattice) delete wArrayLattice; wArrayLattice=0;
822 0 : wArrayLattice = new ArrayLattice<Float>(wGriddedData);
823 0 : wArrayLattice->set(0.0);
824 0 : wLattice=wArrayLattice;
825 :
826 : // Now find the SubLattice corresponding to the image
827 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
828 0 : IPosition stride(4, 1);
829 0 : IPosition trc(blc+image->shape()-stride);
830 0 : LCBox gridBox(blc, trc, gridShape);
831 0 : SubLattice<Complex> gridSub(*arrayLattice, gridBox, true);
832 :
833 : // Do the copy
834 0 : gridSub.copyData(*image);
835 :
836 0 : lattice=arrayLattice;
837 0 : }
838 :
839 0 : AlwaysAssert(lattice, AipsError);
840 0 : AlwaysAssert(wLattice, AipsError);
841 :
842 0 : }
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 0 : void SDGrid::get(VisBuffer& vb, Int row)
1242 : {
1243 0 : LogIO os(LogOrigin("SDGrid", "get"));
1244 :
1245 0 : gridOk(convSupport);
1246 :
1247 : // If row is -1 then we pass through all rows
1248 : Int startRow, endRow, nRow;
1249 0 : if (row==-1) {
1250 0 : nRow=vb.nRow();
1251 0 : startRow=0;
1252 0 : endRow=nRow-1;
1253 0 : 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 0 : 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 0 : if(doConversion_p[vb.spectralWindow()]){
1276 0 : matchChannel(vb.spectralWindow(), vb);
1277 : }
1278 : else{
1279 0 : chanMap.resize();
1280 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
1281 : }
1282 :
1283 : //No point in reading data if its not matching in frequency
1284 0 : if(max(chanMap)==-1)
1285 0 : return;
1286 0 : Cube<Complex> data;
1287 0 : Cube<Int> flags;
1288 0 : getInterpolateArrays(vb, data, flags);
1289 :
1290 : Complex *datStorage;
1291 : Bool isCopy;
1292 0 : datStorage=data.getStorage(isCopy);
1293 : // NOTE: with MS V2.0 the pointing could change per antenna and timeslot
1294 : //
1295 0 : Vector<Int> rowFlags(vb.flagRow().nelements());
1296 0 : rowFlags=0;
1297 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1298 0 : if(vb.flagRow()(rownr)) rowFlags(rownr)=1;
1299 : //single dish yes ?
1300 0 : 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 0 : Matrix<Double> xyPositions(2, endRow-startRow+1);
1352 0 : xyPositions=-1e9;
1353 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1354 0 : if(getXYPos(vb, rownr)) {
1355 0 : xyPositions(0, rownr)=xyPos(0);
1356 0 : xyPositions(1, rownr)=xyPos(1);
1357 : }
1358 : }
1359 :
1360 : Bool del;
1361 : // IPosition s(data.shape());
1362 0 : const IPosition& fs=data.shape();
1363 0 : std::vector<Int> s(fs.begin(), fs.end());
1364 0 : dgridsd(xyPositions.getStorage(del),
1365 : datStorage,
1366 0 : &s[0],
1367 0 : &s[1],
1368 0 : flags.getStorage(del),
1369 0 : rowFlags.getStorage(del),
1370 0 : &s[2],
1371 : &row,
1372 0 : 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 0 : data.putStorage(datStorage, isCopy);
1384 0 : }
1385 0 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1386 0 : }
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()", WHERE);
1471 0 : LogIO logger(logOrigin);
1472 0 : logger << msg.str() << LogIO::SEVERE << LogIO::POST;
1473 :
1474 : // Abort
1475 0 : logger << "Polarization Frame must not change" << LogIO::EXCEPTION;
1476 0 : };
1477 0 : }
1478 :
1479 : } // SDGrid::makeImage helper functions namespace
1480 :
1481 0 : void SDGrid::nextChunk(ROVisibilityIterator &vi) {
1482 : #if defined(SDGRID_PERFS)
1483 0 : StartStop trigger(cNextChunk);
1484 : #endif
1485 0 : vi.nextChunk();
1486 0 : }
1487 :
1488 : // Make a plain straightforward honest-to-FSM image. This returns
1489 : // a complex image, without conversion to Stokes. The representation
1490 : // is that required for the visibilities.
1491 : //----------------------------------------------------------------------
1492 0 : void SDGrid::makeImage(FTMachine::Type inType,
1493 : ROVisibilityIterator& vi,
1494 : ImageInterface<Complex>& theImage,
1495 : Matrix<Float>& weight) {
1496 :
1497 : // Attach visibility buffer (VisBuffer) to visibility iterator (VisibilityIterator)
1498 0 : VisBuffer vb(vi);
1499 :
1500 : // Set the Polarization Representation
1501 : // of image's Stokes Coordinate Axis
1502 : // based on first data in first MS
1503 0 : vi.origin();
1504 0 : const auto imgPolRep = (vb.polFrame() == MSIter::Linear) ?
1505 0 : StokesImageUtil::LINEAR : StokesImageUtil::CIRCULAR;
1506 0 : StokesImageUtil::changeCStokesRep(theImage, imgPolRep);
1507 0 : const auto firstMsName = vb.msName();
1508 :
1509 0 : initializeToSky(theImage, weight, vb);
1510 :
1511 : // Setup SDGrid Cache Manager
1512 0 : const auto onDuty = cacheIsEnabled;
1513 0 : const auto accessMode = cache.isEmpty() ? Cache::AccessMode::WRITE
1514 0 : : Cache::AccessMode::READ;
1515 0 : CacheManager cacheManager(cache, onDuty, accessMode);
1516 :
1517 : // Loop over the visibilities, putting VisBuffers
1518 0 : for (vi.originChunks(); vi.moreChunks(); nextChunk(vi)) {
1519 0 : abortOnPolFrameChange(imgPolRep, firstMsName, vi);
1520 : FTMachine::Type actualType;
1521 : Bool doPSF;
1522 0 : if (vi.newMS()) { // Note: the first MS is a new MS
1523 0 : getParamsForFTMachineType(vi, inType, doPSF, actualType);
1524 0 : handleNewMs(vi, theImage);
1525 : }
1526 0 : for (vi.origin(); vi.more(); vi++) {
1527 0 : setupVisBufferForFTMachineType(actualType, vb);
1528 0 : constexpr Int allVbRows = -1;
1529 0 : put(vb, allVbRows, doPSF, actualType);
1530 : }
1531 : }
1532 0 : finalizeToSky();
1533 :
1534 : // Normalize by dividing out weights, etc.
1535 0 : auto doNormalize = (inType != FTMachine::COVERAGE);
1536 0 : getImage(weight, doNormalize);
1537 :
1538 : // Warning message
1539 0 : if (allEQ(weight, 0.0f)) {
1540 0 : LogIO logger(LogOrigin(name(),"makeImage"));
1541 : logger << LogIO::SEVERE
1542 : << "No useful data in SDGrid: all weights are zero"
1543 0 : << LogIO::POST;
1544 0 : }
1545 0 : }
1546 :
1547 : // Finalize : optionally normalize by weight image
1548 0 : ImageInterface<Complex>& SDGrid::getImage(Matrix<Float>& weights,
1549 : Bool normalize)
1550 : {
1551 0 : AlwaysAssert(lattice, AipsError);
1552 0 : AlwaysAssert(image, AipsError);
1553 :
1554 0 : logIO() << LogOrigin("SDGrid", "getImage") << LogIO::NORMAL;
1555 :
1556 : // execute minmax clipping
1557 0 : clipMinMax();
1558 :
1559 0 : weights.resize(sumWeight.shape());
1560 :
1561 0 : convertArray(weights,sumWeight);
1562 :
1563 : // If the weights are all zero then we cannot normalize
1564 : // otherwise we don't care.
1565 : ///////////////////////
1566 : /*{
1567 : PagedImage<Float> thisScreen(lattice->shape(), image->coordinates(), "TheData");
1568 : LatticeExpr<Float> le(abs(*lattice));
1569 : thisScreen.copyData(le);
1570 : PagedImage<Float> thisScreen2(lattice->shape(), image->coordinates(), "TheWeight");
1571 : LatticeExpr<Float> le2(abs(*wLattice));
1572 : thisScreen2.copyData(le2);
1573 : }*/
1574 : /////////////////////
1575 :
1576 0 : if(normalize) {
1577 0 : if(max(weights)==0.0) {
1578 : //logIO() << LogIO::WARN << "No useful data in SDGrid: weights all zero"
1579 : // << LogIO::POST;
1580 : //image->set(Complex(0.0));
1581 0 : return *image;
1582 : }
1583 : //Timer tim;
1584 : //tim.mark();
1585 : //lattice->copyData((LatticeExpr<Complex>)(iif((*wLattice<=minWeight_p), 0.0,
1586 : // (*lattice)/(*wLattice))));
1587 : //As we are not using disk based lattices anymore the above uses too much memory for
1588 : // ArrayLattice...it does not do a real inplace math but rather mutiple copies
1589 : // seem to be involved thus the less elegant but faster and less memory hog loop
1590 : // below.
1591 :
1592 0 : IPosition pos(4);
1593 0 : for (Int chan=0; chan < nchan; ++chan){
1594 0 : pos[3]=chan;
1595 0 : for( Int pol=0; pol < npol; ++ pol){
1596 0 : pos[2]=pol;
1597 0 : for (Int y=0; y < ny ; ++y){
1598 0 : pos[1]=y;
1599 0 : for (Int x=0; x < nx; ++x){
1600 0 : pos[0]=x;
1601 0 : Float wgt=wGriddedData(pos);
1602 0 : if(wgt > minWeight_p)
1603 0 : griddedData(pos)=griddedData(pos)/wgt;
1604 : else
1605 0 : griddedData(pos)=0.0;
1606 : }
1607 : }
1608 : }
1609 : }
1610 : //tim.show("After normalizing");
1611 0 : }
1612 :
1613 : //if(!isTiled)
1614 : {
1615 : // Now find the SubLattice corresponding to the image
1616 0 : IPosition gridShape(4, nx, ny, npol, nchan);
1617 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
1618 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1619 0 : IPosition stride(4, 1);
1620 0 : IPosition trc(blc+image->shape()-stride);
1621 0 : LCBox gridBox(blc, trc, gridShape);
1622 0 : SubLattice<Complex> gridSub(*arrayLattice, gridBox);
1623 :
1624 : // Do the copy
1625 0 : image->copyData(gridSub);
1626 0 : }
1627 0 : return *image;
1628 : }
1629 :
1630 : // Return weights image
1631 0 : void SDGrid::getWeightImage(ImageInterface<Float>& weightImage, Matrix<Float>& weights)
1632 : {
1633 0 : AlwaysAssert(lattice, AipsError);
1634 0 : AlwaysAssert(image, AipsError);
1635 :
1636 0 : logIO() << LogOrigin("SDGrid", "getWeightImage") << LogIO::NORMAL;
1637 :
1638 0 : weights.resize(sumWeight.shape());
1639 0 : convertArray(weights,sumWeight);
1640 :
1641 0 : weightImage.copyData(*wArrayLattice);
1642 0 : }
1643 :
1644 0 : void SDGrid::ok() {
1645 0 : AlwaysAssert(image, AipsError);
1646 0 : }
1647 :
1648 : // Get the index into the pointing table for this time. Note that the
1649 : // in the pointing table, TIME specifies the beginning of the spanned
1650 : // time range, whereas for the main table, TIME is the centroid.
1651 : // Note that the behavior for multiple matches is not specified! i.e.
1652 : // if there are multiple matches, the index returned depends on the
1653 : // history of previous matches. It is deterministic but not obvious.
1654 : // One could cure this by searching but it would be considerably
1655 : // costlier.
1656 0 : Int SDGrid::getIndex(const MSPointingColumns& mspc, const Double& time,
1657 : const Double& interval, const Int& antid) {
1658 : //Int start=lastIndex_p;
1659 0 : Int start=lastIndexPerAnt_p[antid];
1660 0 : Double tol=interval < 0.0 ? time*DBL_EPSILON : interval/2.0;
1661 : // Search forwards
1662 0 : Int nrows=mspc.time().nrow();
1663 0 : for (Int i=start;i<nrows;i++) {
1664 : // Check for ANTENNA_ID. When antid<0 (default) ANTENNA_ID in POINTING table < 0 is ignored.
1665 : // MS v2 definition indicates ANTENNA_ID in POINING >= 0.
1666 0 : if (antid >= 0 && mspc.antennaId()(i) != antid) {
1667 0 : continue;
1668 : }
1669 :
1670 0 : Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
1671 : // If the interval in the pointing table is negative, use the last
1672 : // entry. Note that this may be invalid (-1) but in that case
1673 : // the calling routine will generate an error
1674 0 : Double mspc_interval = mspc.interval()(i);
1675 :
1676 0 : if(mspc_interval<0.0) {
1677 : //return lastIndex_p;
1678 0 : return lastIndexPerAnt_p[antid];
1679 : }
1680 : // Pointing table interval is specified so we have to do a match
1681 : else {
1682 : // Is the midpoint of this pointing table entry within the specified
1683 : // tolerance of the main table entry?
1684 0 : if(abs(midpoint-time) <= (mspc_interval/2.0+tol)) {
1685 : //lastIndex_p=i;
1686 0 : lastIndexPerAnt_p[antid]=i;
1687 0 : return i;
1688 : }
1689 : }
1690 : }
1691 : // Search backwards
1692 0 : for (Int i=start;i>=0;i--) {
1693 0 : if (antid >= 0 && mspc.antennaId()(i) != antid) {
1694 0 : continue;
1695 : }
1696 0 : Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
1697 0 : Double mspc_interval = mspc.interval()(i);
1698 0 : if(mspc_interval<0.0) {
1699 : //return lastIndex_p;
1700 0 : return lastIndexPerAnt_p[antid];
1701 : }
1702 : // Pointing table interval is specified so we have to do a match
1703 : else {
1704 : // Is the midpoint of this pointing table entry within the specified
1705 : // tolerance of the main table entry?
1706 0 : if(abs(midpoint-time) <= (mspc_interval/2.0+tol)) {
1707 : //lastIndex_p=i;
1708 0 : lastIndexPerAnt_p[antid]=i;
1709 0 : return i;
1710 : }
1711 : }
1712 : }
1713 : // No match!
1714 0 : return -1;
1715 : }
1716 :
1717 0 : Bool SDGrid::getXYPos(const VisBuffer& vb, Int row) {
1718 :
1719 : // Cache control
1720 0 : if (cacheIsEnabled and cache.isReadable()) {
1721 0 : cache.loadRowPixel();
1722 0 : return rowPixel.isValid;
1723 : }
1724 0 : const auto onDuty = cacheIsEnabled and cache.isWriteable();
1725 0 : CacheWriter cacheWriter(cache, onDuty);
1726 :
1727 : // Until we manage to compute a valid one ...
1728 0 : rowPixel.isValid = false;
1729 :
1730 : // Select the POINTING table (columns) we'll work with.
1731 0 : const auto haveConvertedColumn = ramPointingTable.nrow() > 0;
1732 0 : const MSPointingColumns& pointingColumns = haveConvertedColumn ?
1733 0 : *ramPointingColumnsPtr
1734 0 : : vb.msColumns().pointing();
1735 :
1736 0 : const auto nPointings = pointingColumns.nrow();
1737 0 : const auto havePointings = (nPointings >= 1);
1738 :
1739 : // We'll need to call these many times, so let's call them once for good
1740 0 : const auto rowTime = vb.time()(row);
1741 0 : const auto rowTimeInterval = vb.timeInterval()(row);
1742 0 : const auto rowAntenna1 = vb.antenna1()(row);
1743 :
1744 : // 1. Try to find the index of a pointing recorded:
1745 : // - for the antenna of specified row,
1746 : // - at a time close enough to the time at which
1747 : // data of specified row was taken using that antenna
1748 0 : constexpr Int invalidIndex = -1;
1749 0 : Int pointingIndex = invalidIndex;
1750 0 : if (havePointings) {
1751 : #if defined(SDGRID_PERFS)
1752 0 : StartStop trigger(cSearchValidPointing);
1753 : #endif
1754 : // if (vb.newMS()) vb.newMS does not work well using msid
1755 : // Note about above comment:
1756 : // - vb.newMS probably works well
1757 : // - but if the calling code is iterating over the rows of a subchunk
1758 : // vb.newMS returns true for all rows belonging to the first subchunk
1759 : // of the first chunk of a new MS.
1760 :
1761 : // ???
1762 : // What if vb changed since we were last called ?
1763 : // What if the calling code calls put and get, with different VisBuffers ?
1764 0 : if (vb.msId() != msId_p) {
1765 0 : lastIndex_p = 0; // no longer used
1766 0 : if (lastIndexPerAnt_p.nelements() < (size_t)vb.numberAnt()) {
1767 0 : lastIndexPerAnt_p.resize(vb.numberAnt());
1768 : }
1769 0 : lastIndexPerAnt_p = 0;
1770 0 : msId_p = vb.msId();
1771 0 : lastAntID_p = -1;
1772 : }
1773 :
1774 : // Try to locate a pointing verifying for specified row:
1775 : // | POINTING.TIME - MAIN.TIME | <= 0.5*(POINTING.INTERVAL + tolerance)
1776 :
1777 : // Try first using a tiny tolerance
1778 0 : constexpr Double useTinyTolerance = -1.0;
1779 0 : pointingIndex = getIndex(pointingColumns, rowTime, useTinyTolerance , rowAntenna1);
1780 :
1781 0 : auto foundPointing = (pointingIndex >= 0);
1782 0 : if (not foundPointing) {
1783 : // Try again using tolerance = MAIN.INTERVAL
1784 0 : pointingIndex = getIndex(pointingColumns, rowTime, rowTimeInterval, rowAntenna1);
1785 0 : foundPointing = (pointingIndex >= 0);
1786 : }
1787 :
1788 : // Making the implicit type conversion explicit.
1789 : // Conversion is safe because it occurs only when pointingIndex >= 0.
1790 0 : const auto foundValidPointing = (foundPointing and (static_cast<rownr_t>(pointingIndex) < nPointings));
1791 :
1792 0 : if (not foundValidPointing) {
1793 0 : ostringstream o;
1794 0 : o << "Failed to find pointing information for time "
1795 0 : << MVTime(rowTime/86400.0) << " : Omitting this point";
1796 :
1797 0 : logIO_p << LogIO::DEBUGGING << String(o) << LogIO::POST;
1798 :
1799 0 : return rowPixel.isValid;
1800 0 : }
1801 0 : }
1802 :
1803 : // 2. At this stage we have:
1804 : // * either no pointings and an invalid pointingIndex
1805 : // * or pointings and a valid pointingIndex.
1806 : // Decide now if we need to interpolate antenna's pointing direction
1807 : // at data-taking time:
1808 : // we'll do so when data is sampled faster than pointings are recorded
1809 0 : Bool needInterpolation = False;
1810 0 : if (havePointings) {
1811 0 : const auto pointingInterval = pointingColumns.interval()(pointingIndex);
1812 0 : if (rowTimeInterval < pointingInterval) needInterpolation = True;
1813 : }
1814 0 : const auto mustInterpolate = havePointings && needInterpolation;
1815 :
1816 : // 3. Create interpolator if needed
1817 0 : if (mustInterpolate) {
1818 0 : if (not isSplineInterpolationReady) {
1819 : #if defined(SDGRID_PERFS)
1820 0 : StartStop trigger(cComputeSplines);
1821 : #endif
1822 : const auto nAntennas = static_cast<size_t>(
1823 0 : vb.msColumns().antenna().nrow()
1824 : );
1825 0 : interpolator = new SDPosInterpolator(
1826 : pointingColumns,
1827 0 : pointingDirCol_p,
1828 : nAntennas
1829 0 : );
1830 0 : isSplineInterpolationReady = true;
1831 0 : } else {
1832 : // We have an interpolator. Re-use it if possible.
1833 0 : const auto canReuseInterpolator = interpolator->inTimeRange(rowTime, rowAntenna1);
1834 0 : if (not canReuseInterpolator) {
1835 : // setup spline interpolator for the current dataset (CAS-11261, 2018/5/22 WK)
1836 : // delete and re-create it
1837 0 : delete interpolator;
1838 0 : interpolator = 0;
1839 : const auto nAntennas = static_cast<size_t>(
1840 0 : vb.msColumns().antenna().nrow()
1841 : );
1842 0 : interpolator = new SDPosInterpolator(
1843 : pointingColumns,
1844 0 : pointingDirCol_p,
1845 : nAntennas
1846 0 : );
1847 : }
1848 : }
1849 : }
1850 :
1851 : // 4. Create the direction conversion machine if needed
1852 :
1853 0 : if ( pointingDirCol_p == "SOURCE_OFFSET" or
1854 0 : pointingDirCol_p == "POINTING_OFFSET" ) {
1855 : // It makes no sense to track in offset coordinates...
1856 : // hopefully the user sets the image coords right
1857 0 : fixMovingSource_p = false;
1858 : }
1859 :
1860 0 : const auto needDirectionConverter = (
1861 0 : not havePointings or not haveConvertedColumn or fixMovingSource_p
1862 : );
1863 :
1864 0 : if (needDirectionConverter) {
1865 0 : if (not pointingToImage) {
1866 : // Set the frame
1867 : const auto & rowAntenna1Position =
1868 0 : vb.msColumns().antenna().positionMeas()(rowAntenna1);
1869 : // set dummy time stamp 1 day before rowTime
1870 0 : const MEpoch dummyEpoch(Quantity(rowTime - 86400.0, "s"));
1871 0 : mFrame_p = MeasFrame(dummyEpoch, rowAntenna1Position);
1872 :
1873 : // Remember antenna id for next call,
1874 : // which may be done using a different VisBuffer ...
1875 0 : lastAntID_p = rowAntenna1;
1876 :
1877 : // Compute the "model" required to setup the conversion machine
1878 0 : if (havePointings) {
1879 0 : worldPosMeas = mustInterpolate ?
1880 : directionMeas(pointingColumns, pointingIndex, rowTime)
1881 0 : : directionMeas(pointingColumns, pointingIndex);
1882 : } else {
1883 : // Without pointings, this sets the direction to the phase center
1884 0 : worldPosMeas = vb.direction1()(row);
1885 : }
1886 :
1887 : // Make a machine to convert from the worldPosMeas to the output
1888 : // Direction Measure type for the relevant frame
1889 0 : MDirection::Ref outRef(directionCoord.directionType(), mFrame_p);
1890 0 : pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
1891 0 : if (not pointingToImage) {
1892 0 : logIO_p << "Cannot make direction conversion machine" << LogIO::EXCEPTION;
1893 : }
1894 : // Perform 1 dummy direction conversion to clear values
1895 : // cached in static variables of casacore functions like MeasTable::dUT1
1896 0 : MDirection _dir_tmp = (*pointingToImage)();
1897 0 : }
1898 : }
1899 :
1900 : // 5. Update the frame holding the measurements for this row
1901 0 : const MEpoch rowEpoch(Quantity(rowTime, "s"));
1902 : // ---- Always reset the epoch
1903 : {
1904 : #if defined(SDGRID_PERFS)
1905 0 : StartStop trigger(cResetFrame);
1906 : #endif
1907 0 : mFrame_p.resetEpoch(rowEpoch);
1908 0 : }
1909 : // ---- Reset antenna position only if antenna changed since we were last called
1910 0 : if (lastAntID_p != rowAntenna1) {
1911 : // Debug messages
1912 0 : if (lastAntID_p == -1) {
1913 : // antenna ID is unset
1914 0 : logIO_p << LogIO::DEBUGGING
1915 : << "updating antenna position for conversion: new MS ID " << msId_p
1916 0 : << ", antenna ID " << rowAntenna1 << LogIO::POST;
1917 : } else {
1918 0 : logIO_p << LogIO::DEBUGGING
1919 : << "updating antenna position for conversion: MS ID " << msId_p
1920 : << ", last antenna ID " << lastAntID_p
1921 0 : << ", new antenna ID " << rowAntenna1 << LogIO::POST;
1922 : }
1923 : const auto & rowAntenna1Position =
1924 0 : vb.msColumns().antenna().positionMeas()(rowAntenna1);
1925 : {
1926 : #if defined(SDGRID_PERFS)
1927 0 : StartStop trigger(cResetFrame);
1928 : #endif
1929 0 : mFrame_p.resetPosition(rowAntenna1Position);
1930 0 : }
1931 : // Remember antenna id for next call,
1932 : // which may be done using a different VisBuffer ...
1933 0 : lastAntID_p = rowAntenna1;
1934 0 : }
1935 :
1936 : // 6. Compute user-specified column direction at data-taking time,
1937 : // in image's direction reference frame
1938 0 : if (havePointings) {
1939 0 : if (mustInterpolate) {
1940 : #if defined(SDGRID_PERFS)
1941 0 : cInterpolateDirection.start();
1942 : #endif
1943 : const auto interpolatedDirection =
1944 0 : directionMeas(pointingColumns, pointingIndex, rowTime);
1945 : #if defined(SDGRID_PERFS)
1946 0 : cInterpolateDirection.stop();
1947 0 : cConvertDirection.start();
1948 : #endif
1949 : worldPosMeas = haveConvertedColumn ?
1950 : interpolatedDirection
1951 0 : : (*pointingToImage)(interpolatedDirection);
1952 : #if defined(SDGRID_PERFS)
1953 0 : cConvertDirection.stop();
1954 : #endif
1955 :
1956 : // Debug stuff
1957 : //Vector<Double> newdirv = newdir.getAngle("rad").getValue();
1958 : //cerr<<"dir0="<<newdirv(0)<<endl;
1959 :
1960 : //fprintf(pfile,"%.8f %.8f \n", newdirv(0), newdirv(1));
1961 : //printf("%lf %lf \n", newdirv(0), newdirv(1));
1962 0 : } else {
1963 0 : const auto columnDirection = directionMeas(pointingColumns, pointingIndex);
1964 : worldPosMeas = haveConvertedColumn ?
1965 : columnDirection
1966 0 : : (*pointingToImage)(columnDirection);
1967 0 : }
1968 : } else {
1969 : // Without pointings, this converts the direction of the phase center
1970 0 : worldPosMeas = (*pointingToImage)(vb.direction1()(row));
1971 : }
1972 :
1973 : // 7. Convert world position coordinates to image pixel coordinates
1974 : {
1975 : #if defined(SDGRID_PERFS)
1976 0 : StartStop trigger(cComputeDirectionPixel);
1977 : #endif
1978 :
1979 0 : rowPixel.isValid = directionCoord.toPixel(xyPos, worldPosMeas);
1980 0 : }
1981 :
1982 0 : if (not rowPixel.isValid) {
1983 0 : logIO_p << "Failed to find a pixel for pointing direction of "
1984 0 : << MVTime(worldPosMeas.getValue().getLong("rad")).string(MVTime::TIME)
1985 0 : << ", " << MVAngle(worldPosMeas.getValue().getLat("rad")).string(MVAngle::ANGLE)
1986 0 : << LogIO::WARN << LogIO::POST;
1987 0 : return rowPixel.isValid;
1988 : }
1989 :
1990 : // 8. Handle moving sources
1991 0 : if (fixMovingSource_p) {
1992 : #if defined(SDGRID_PERFS)
1993 0 : StartStop trigger(cHandleMovingSource);
1994 : #endif
1995 0 : if (xyPosMovingOrig_p.nelements() < 2) {
1996 0 : directionCoord.toPixel(xyPosMovingOrig_p, firstMovingDir_p);
1997 : }
1998 : //via HADEC or AZEL for parallax of near sources
1999 0 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
2000 0 : MDirection tmphadec = MDirection::Convert(movingDir_p, outref1)();
2001 0 : MDirection actSourceDir = (*pointingToImage)(tmphadec);
2002 0 : Vector<Double> actPix;
2003 0 : directionCoord.toPixel(actPix, actSourceDir);
2004 :
2005 : //cout << row << " scan " << vb.scan()(row) << "xyPos " << xyPos
2006 : // << " xyposmovorig " << xyPosMovingOrig_p << " actPix " << actPix << endl;
2007 :
2008 0 : xyPos = xyPos + xyPosMovingOrig_p - actPix;
2009 0 : }
2010 :
2011 0 : return rowPixel.isValid;
2012 0 : }
2013 :
2014 0 : MDirection SDGrid::directionMeas(const MSPointingColumns& mspc, const Int& index){
2015 0 : if (pointingDirCol_p == "TARGET") {
2016 0 : return mspc.targetMeas(index);
2017 0 : } else if (pointingDirCol_p == "POINTING_OFFSET") {
2018 0 : if (!mspc.pointingOffsetMeasCol().isNull()) {
2019 0 : return mspc.pointingOffsetMeas(index);
2020 : }
2021 0 : cerr << "No PONTING_OFFSET column in POINTING table" << endl;
2022 0 : } else if (pointingDirCol_p == "SOURCE_OFFSET") {
2023 0 : if (!mspc.sourceOffsetMeasCol().isNull()) {
2024 0 : return mspc.sourceOffsetMeas(index);
2025 : }
2026 0 : cerr << "No SOURCE_OFFSET column in POINTING table" << endl;
2027 0 : } else if (pointingDirCol_p == "ENCODER") {
2028 0 : if (!mspc.encoderMeas().isNull()) {
2029 0 : return mspc.encoderMeas()(index);
2030 : }
2031 0 : cerr << "No ENCODER column in POINTING table" << endl;
2032 : }
2033 :
2034 : //default return this
2035 0 : return mspc.directionMeas(index);
2036 : }
2037 :
2038 : // for the cases, interpolation of the pointing direction requires
2039 : // when data sampling rate higher than the pointing data recording
2040 : // (e.g. fast OTF)
2041 0 : MDirection SDGrid::directionMeas(const MSPointingColumns& mspc, const Int& index,
2042 : const Double& time){
2043 : //spline interpolation
2044 0 : if (isSplineInterpolationReady) {
2045 0 : Int antid = mspc.antennaId()(index);
2046 0 : if (interpolator->doSplineInterpolation(antid)) {
2047 0 : return interpolator->interpolateDirectionMeasSpline(mspc, time, index, antid);
2048 : }
2049 : }
2050 :
2051 : //linear interpolation (as done before CAS-7911)
2052 : Int index1, index2;
2053 0 : if (time < mspc.time()(index)) {
2054 0 : if (index > 0) {
2055 0 : index1 = index-1;
2056 0 : index2 = index;
2057 0 : } else if (index == 0) {
2058 0 : index1 = index;
2059 0 : index2 = index+1;
2060 : }
2061 : } else {
2062 0 : if (index < Int(mspc.nrow()-1)) {
2063 0 : index1 = index;
2064 0 : index2 = index+1;
2065 0 : } else if (index == Int(mspc.nrow()-1) || (mspc.time()(index)-mspc.time()(index+1))>2*mspc.interval()(index)) {
2066 0 : index1 = index-1;
2067 0 : index2 = index;
2068 : }
2069 : }
2070 0 : return interpolateDirectionMeas(mspc, time, index, index1, index2);
2071 : }
2072 :
2073 0 : MDirection SDGrid::interpolateDirectionMeas(const MSPointingColumns& mspc,
2074 : const Double& time,
2075 : const Int& index,
2076 : const Int& indx1,
2077 : const Int& indx2){
2078 0 : Vector<Double> dir1,dir2;
2079 0 : Vector<Double> newdir(2),scanRate(2);
2080 : Double dLon, dLat;
2081 : Double ftime,ftime2,ftime1,dtime;
2082 0 : MDirection newDirMeas;
2083 0 : MDirection::Ref rf;
2084 : Bool isfirstRefPt;
2085 :
2086 0 : if (indx1 == index) {
2087 0 : isfirstRefPt = true;
2088 : } else {
2089 0 : isfirstRefPt = false;
2090 : }
2091 :
2092 0 : if (pointingDirCol_p == "TARGET") {
2093 0 : dir1 = mspc.targetMeas(indx1).getAngle("rad").getValue();
2094 0 : dir2 = mspc.targetMeas(indx2).getAngle("rad").getValue();
2095 0 : } else if (pointingDirCol_p == "POINTING_OFFSET") {
2096 0 : if (!mspc.pointingOffsetMeasCol().isNull()) {
2097 0 : dir1 = mspc.pointingOffsetMeas(indx1).getAngle("rad").getValue();
2098 0 : dir2 = mspc.pointingOffsetMeas(indx2).getAngle("rad").getValue();
2099 : } else {
2100 0 : cerr << "No PONTING_OFFSET column in POINTING table" << endl;
2101 : }
2102 0 : } else if (pointingDirCol_p == "SOURCE_OFFSET") {
2103 0 : if (!mspc.sourceOffsetMeasCol().isNull()) {
2104 0 : dir1 = mspc.sourceOffsetMeas(indx1).getAngle("rad").getValue();
2105 0 : dir2 = mspc.sourceOffsetMeas(indx2).getAngle("rad").getValue();
2106 : } else {
2107 0 : cerr << "No SOURCE_OFFSET column in POINTING table" << endl;
2108 : }
2109 0 : } else if (pointingDirCol_p == "ENCODER") {
2110 0 : if (!mspc.encoderMeas().isNull()) {
2111 0 : dir1 = mspc.encoderMeas()(indx1).getAngle("rad").getValue();
2112 0 : dir2 = mspc.encoderMeas()(indx2).getAngle("rad").getValue();
2113 : } else {
2114 0 : cerr << "No ENCODER column in POINTING table" << endl;
2115 : }
2116 : } else {
2117 0 : dir1 = mspc.directionMeas(indx1).getAngle("rad").getValue();
2118 0 : dir2 = mspc.directionMeas(indx2).getAngle("rad").getValue();
2119 : }
2120 :
2121 0 : dLon = dir2(0) - dir1(0);
2122 0 : dLat = dir2(1) - dir1(1);
2123 0 : ftime = floor(mspc.time()(indx1));
2124 0 : ftime2 = mspc.time()(indx2) - ftime;
2125 0 : ftime1 = mspc.time()(indx1) - ftime;
2126 0 : dtime = ftime2 - ftime1;
2127 0 : scanRate(0) = dLon/dtime;
2128 0 : scanRate(1) = dLat/dtime;
2129 : //scanRate(0) = dir2(0)/dtime-dir1(0)/dtime;
2130 : //scanRate(1) = dir2(1)/dtime-dir1(1)/dtime;
2131 : //Double delT = mspc.time()(index)-time;
2132 : //cerr<<"index="<<index<<" dLat="<<dLat<<" dtime="<<dtime<<" delT="<< delT<<endl;
2133 : //cerr<<"deldirlat="<<scanRate(1)*fabs(delT)<<endl;
2134 0 : if (isfirstRefPt) {
2135 0 : newdir(0) = dir1(0)+scanRate(0)*fabs(mspc.time()(index)-time);
2136 0 : newdir(1) = dir1(1)+scanRate(1)*fabs(mspc.time()(index)-time);
2137 0 : rf = mspc.directionMeas(indx1).getRef();
2138 : } else {
2139 0 : newdir(0) = dir2(0)-scanRate(0)*fabs(mspc.time()(index)-time);
2140 0 : newdir(1) = dir2(1)-scanRate(1)*fabs(mspc.time()(index)-time);
2141 0 : rf = mspc.directionMeas(indx2).getRef();
2142 : }
2143 : //default return this
2144 0 : Quantity rDirLon(newdir(0), "rad");
2145 0 : Quantity rDirLat(newdir(1), "rad");
2146 0 : newDirMeas = MDirection(rDirLon, rDirLat, rf);
2147 : //cerr<<"newDirMeas rf="<<newDirMeas.getRefString()<<endl;
2148 : //return mspc.directionMeas(index);
2149 0 : return newDirMeas;
2150 0 : }
2151 :
2152 0 : void SDGrid::pickWeights(const VisBuffer& vb, Matrix<Float>& weight){
2153 : #if defined(SDGRID_PERFS)
2154 0 : StartStop trigger(cPickWeights);
2155 : #endif
2156 : //break reference
2157 0 : weight.resize();
2158 :
2159 0 : if (useImagingWeight_p) {
2160 0 : weight.reference(vb.imagingWeight());
2161 : } else {
2162 0 : const Cube<Float> weightSpec(vb.weightSpectrum());
2163 0 : weight.resize(vb.nChannel(), vb.nRow());
2164 :
2165 : // CAS-9957 correct weight propagation from linear/circular correlations to Stokes I
2166 0 : const auto toStokesWeight = [](float weight0, float weight1) {
2167 0 : const auto denominator = weight0 + weight1;
2168 0 : const auto numerator = weight0 * weight1;
2169 0 : constexpr float fmin = std::numeric_limits<float>::min();
2170 0 : return abs(denominator) < fmin ? 0.0f : 4.0f * numerator / denominator;
2171 : };
2172 :
2173 0 : if (weightSpec.nelements() == 0) {
2174 0 : const auto &weightMat = vb.weightMat();
2175 0 : const ssize_t npol = weightMat.shape()(0);
2176 0 : if (npol == 1) {
2177 0 : const auto weight0 = weightMat.row(0);
2178 0 : for (int k = 0; k < vb.nRow(); ++k) {
2179 0 : weight.column(k).set(weight0(k));
2180 : }
2181 0 : } else if (npol == 2) {
2182 0 : const auto weight0 = weightMat.row(0);
2183 0 : const auto weight1 = weightMat.row(1);
2184 0 : for (int k = 0; k < vb.nRow(); ++k) {
2185 : //cerr << "nrow " << vb.nRow() << " " << weight.shape() << " " << weight.column(k).shape() << endl;
2186 0 : weight.column(k).set(toStokesWeight(weight0(k), weight1(k)));
2187 : }
2188 0 : } else {
2189 : // It seems current code doesn't support 4 pol case. So, give up
2190 : // processing such data to avoid producing unintended result
2191 0 : throw AipsError("Imaging full-Stokes data (npol=4) is not supported.");
2192 : }
2193 : } else {
2194 0 : const ssize_t npol = weightSpec.shape()(0);
2195 0 : if (npol == 1) {
2196 0 : weight = weightSpec.yzPlane(0);
2197 0 : } else if (npol == 2) {
2198 0 : const auto weight0 = weightSpec.yzPlane(0);
2199 0 : const auto weight1 = weightSpec.yzPlane(1);
2200 0 : for (int k = 0; k < vb.nRow(); ++k) {
2201 0 : for (int chan = 0; chan < vb.nChannel(); ++chan) {
2202 0 : weight(chan, k) = toStokesWeight(weight0(chan, k), weight1(chan, k));
2203 : }
2204 : }
2205 0 : } else {
2206 : // It seems current code doesn't support 4 pol case. So, give up
2207 : // processing such data to avoid producing unintended result
2208 0 : throw AipsError("Imaging full-Stokes data (npol=4) is not supported.");
2209 : }
2210 : }
2211 0 : }
2212 0 : }
2213 :
2214 0 : void SDGrid::clipMinMax() {
2215 0 : if (clipminmax_) {
2216 : Bool gmin_b, gmax_b, wmin_b, wmax_b, np_b;
2217 0 : const auto *gmin_p = gmin_.getStorage(gmin_b);
2218 0 : const auto *gmax_p = gmax_.getStorage(gmax_b);
2219 0 : const auto *wmin_p = wmin_.getStorage(wmin_b);
2220 0 : const auto *wmax_p = wmax_.getStorage(wmax_b);
2221 0 : const auto *np_p = npoints_.getStorage(np_b);
2222 :
2223 : Bool data_b, weight_b, sumw_b;
2224 0 : auto data_p = griddedData.getStorage(data_b);
2225 0 : auto weight_p = wGriddedData.getStorage(weight_b);
2226 0 : auto sumw_p = sumWeight.getStorage(sumw_b);
2227 :
2228 0 : auto arrayShape = griddedData.shape();
2229 0 : size_t num_xy = arrayShape.getFirst(2).product();
2230 0 : size_t num_polchan = arrayShape.getLast(2).product();
2231 0 : for (size_t i = 0; i < num_xy; ++i) {
2232 0 : for (size_t j = 0; j < num_polchan; ++j) {
2233 0 : auto k = i * num_polchan + j;
2234 0 : if (np_p[k] == 1) {
2235 0 : auto wt = wmin_p[k];
2236 0 : data_p[k] = wt * gmin_p[k];
2237 0 : weight_p[k] = wt;
2238 0 : sumw_p[j] += wt;
2239 0 : } else if (np_p[k] == 2) {
2240 0 : auto wt = wmin_p[k];
2241 0 : data_p[k] = wt * gmin_p[k];
2242 0 : weight_p[k] = wt;
2243 0 : sumw_p[j] += wt;
2244 0 : wt = wmax_p[k];
2245 0 : data_p[k] += wt * gmax_p[k];
2246 0 : weight_p[k] += wt;
2247 0 : sumw_p[j] += wt;
2248 : }
2249 : }
2250 : }
2251 :
2252 0 : wGriddedData.putStorage(weight_p, weight_b);
2253 0 : griddedData.putStorage(data_p, data_b);
2254 0 : sumWeight.putStorage(sumw_p, sumw_b);
2255 :
2256 0 : npoints_.freeStorage(np_p, np_b);
2257 0 : wmax_.freeStorage(wmax_p, wmax_b);
2258 0 : wmin_.freeStorage(wmin_p, wmin_b);
2259 0 : gmax_.freeStorage(gmax_p, gmax_b);
2260 0 : gmin_.freeStorage(gmin_p, gmin_b);
2261 0 : }
2262 0 : }
2263 :
2264 0 : SDGrid::MaskedPixelRef::MaskedPixelRef(Vector<Double>& xyIn, Bool isValidIn)
2265 0 : : xy {xyIn},
2266 0 : isValid {isValidIn}
2267 0 : {}
2268 :
2269 : SDGrid::MaskedPixelRef&
2270 0 : SDGrid::MaskedPixelRef::operator=(const SDGrid::MaskedPixelRef &other) {
2271 0 : xy = other.xy;
2272 0 : isValid = other.isValid;
2273 0 : return *this;
2274 : }
2275 :
2276 0 : SDGrid::MaskedPixel::MaskedPixel(Double xIn, Double yIn, Bool isValidIn)
2277 0 : : x {xIn},
2278 0 : y {yIn},
2279 0 : isValid {isValidIn}
2280 0 : {}
2281 :
2282 0 : SDGrid::Cache::Cache(SDGrid &parent)
2283 0 : : sdgrid {parent},
2284 0 : isOpened {false},
2285 0 : accessMode {AccessMode::READ},
2286 0 : canRead {false},
2287 0 : canWrite {false},
2288 0 : inputPixel {sdgrid.rowPixel},
2289 0 : outputPixel {sdgrid.rowPixel}
2290 0 : {}
2291 :
2292 0 : const casacore::String& SDGrid::Cache::className() {
2293 0 : static casacore::String className_ = "SDGrid::Cache";
2294 0 : return className_;
2295 : }
2296 :
2297 0 : SDGrid::Cache& SDGrid::Cache::operator=(const Cache &other) {
2298 0 : sdgrid = other.sdgrid;
2299 0 : msCaches = other.msCaches;
2300 0 : isOpened = false;
2301 0 : canRead = false;
2302 0 : canWrite = false;
2303 : // inputPixel = sdgrid.rowPixel;
2304 0 : inputPixel.xy = sdgrid.rowPixel.xy;
2305 : //inputPixel.isValid = sdgrid.rowPixel.isValid;
2306 : // <=>
2307 0 : msPixels = nullptr;
2308 0 : outputPixel = other.outputPixel;
2309 0 : msCacheReadIterator = MsCaches::const_iterator();
2310 0 : pixelReadIterator = Pixels::const_iterator();
2311 0 : return *this;
2312 : }
2313 :
2314 : void
2315 0 : SDGrid::Cache::open(AccessMode accessModeIn) {
2316 0 : if (isOpened) {
2317 0 : LogIO logger(LogOrigin(className(), "open", WHERE));
2318 0 : logger << "BUG: Opened " << className() << " was re-opened before being closed." << LogIO::EXCEPTION;
2319 0 : }
2320 0 : isOpened = true;
2321 0 : accessMode = accessModeIn;
2322 0 : canRead = accessMode == AccessMode::READ;
2323 0 : canWrite = accessMode == AccessMode::WRITE;
2324 0 : if (outputPixel.xy.size() < 2) {
2325 0 : outputPixel.xy.resize(2);
2326 0 : outputPixel.isValid = false;
2327 : }
2328 0 : rewind();
2329 0 : }
2330 :
2331 : void
2332 0 : SDGrid::Cache::rewind() {
2333 : // Writer
2334 0 : msPixels = nullptr;
2335 :
2336 : // Reader
2337 0 : msCacheReadIterator = msCaches.cbegin();
2338 0 : pixelReadIterator = Pixels::const_iterator();
2339 0 : }
2340 :
2341 : void
2342 0 : SDGrid::Cache::close() {
2343 0 : if (not isOpened) {
2344 0 : LogIO logger(LogOrigin(className(), "close", WHERE));
2345 0 : logger << "BUG: Closed " << className() << " was re-closed before being opened." << LogIO::EXCEPTION;
2346 0 : }
2347 0 : if (isWriteable()) {
2348 : // Make sure we have 1 pixel per row
2349 0 : for (const auto & msCache : msCaches) {
2350 0 : if (not msCache.isConsistent()) {
2351 0 : LogIO logger(LogOrigin(className(), "close", WHERE));
2352 0 : const auto didOverflow = msCache.pixels.size() > msCache.nRows;
2353 0 : const auto overOrUnder = didOverflow ? "over" : "under";
2354 : logger << "BUG: Cache " << overOrUnder << "flow:"
2355 0 : << " nRows=" << msCache.nRows
2356 : << " != nPixels=" << msCache.pixels.size()
2357 0 : << " for: " << msCache.msPath
2358 0 : << " with selection: " << msCache.msTableName
2359 0 : << LogIO::EXCEPTION;
2360 0 : }
2361 : }
2362 : }
2363 : // Reset state
2364 0 : isOpened = false;
2365 0 : accessMode = AccessMode::READ;
2366 0 : canRead = false;
2367 0 : canWrite = false;
2368 0 : }
2369 :
2370 : Bool
2371 0 : SDGrid::Cache::isReadable() const {
2372 0 : return isOpened and canRead;
2373 : }
2374 :
2375 : Bool
2376 0 : SDGrid::Cache::isWriteable() const {
2377 0 : return isOpened and canWrite;
2378 : }
2379 :
2380 : Bool
2381 0 : SDGrid::Cache::isEmpty() const {
2382 0 : return msCaches.empty();
2383 : }
2384 :
2385 : void
2386 0 : SDGrid::Cache::clear() {
2387 0 : msCaches.clear();
2388 0 : }
2389 :
2390 0 : SDGrid::Cache::MsCache::MsCache(const String& msPathIn, const String& msTableNameIn, rownr_t nRowsIn)
2391 0 : : msPath {msPathIn},
2392 0 : msTableName {msTableNameIn},
2393 0 : nRows {nRowsIn}
2394 : {
2395 0 : pixels.reserve(nRows);
2396 0 : }
2397 :
2398 0 : Bool SDGrid::Cache::MsCache::isConsistent() const {
2399 0 : return pixels.size() == nRows;
2400 : }
2401 :
2402 : void
2403 0 : SDGrid::Cache::newMS(const MeasurementSet& ms) {
2404 0 : LogIO os(LogOrigin(className(),"newMS"));
2405 0 : const auto msPath = ms.getPartNames()[0];
2406 :
2407 0 : if (isWriteable()) {
2408 0 : os << "Will compute and cache spectra pixels coordinates for: " << msPath << LogIO::POST;
2409 0 : msCaches.emplace_back(msPath, ms.tableName(), ms.nrow());
2410 0 : msPixels = &(msCaches.back().pixels);
2411 0 : return;
2412 : }
2413 :
2414 0 : if (isReadable()) {
2415 0 : os << "Will load cached spectra pixels coordinates for: " << msPath << LogIO::POST;
2416 0 : if (msCacheReadIterator == msCaches.cend()) {
2417 0 : os << "BUG! Cached data missing for: " << msPath << LogIO::EXCEPTION;
2418 : }
2419 0 : const auto & pixelsToLoad = msCacheReadIterator->pixels;
2420 0 : if (pixelsToLoad.size() != ms.nrow()) {
2421 : os << "BUG! Cached data size mismatch for: " << msPath
2422 : << " : nRows: " << ms.nrow()
2423 0 : << " nCachedRows: " << pixelsToLoad.size() << LogIO::EXCEPTION;
2424 : }
2425 0 : if ( msCacheReadIterator->msTableName != ms.tableName()) {
2426 : os << "BUG! Cached data was probably for a different selection of: " << msPath
2427 : << " current selection: " << ms.tableName()
2428 0 : << " cached selection: " << msCacheReadIterator->msTableName << LogIO::EXCEPTION;
2429 : }
2430 0 : pixelReadIterator = pixelsToLoad.cbegin();
2431 0 : ++msCacheReadIterator;
2432 : }
2433 0 : }
2434 :
2435 : void
2436 0 : SDGrid::Cache::storeRowPixel() {
2437 0 : msPixels->emplace_back(
2438 0 : inputPixel.xy[0],
2439 0 : inputPixel.xy[1],
2440 0 : inputPixel.isValid
2441 : );
2442 0 : }
2443 :
2444 : void
2445 0 : SDGrid::Cache::loadRowPixel() {
2446 0 : const auto & pixel = *pixelReadIterator;
2447 0 : outputPixel.xy[0] = pixel.x;
2448 0 : outputPixel.xy[1] = pixel.y;
2449 0 : outputPixel.isValid = pixel.isValid;
2450 0 : ++pixelReadIterator;
2451 0 : }
2452 :
2453 0 : SDGrid::CacheManager::CacheManager(Cache& cacheIn,
2454 0 : Bool onDutyIn, Cache::AccessMode accessModeIn)
2455 0 : : cache {cacheIn},
2456 0 : onDuty {onDutyIn},
2457 0 : accessMode {accessModeIn}
2458 : {
2459 0 : if (onDuty) {
2460 0 : cache.open(accessMode);
2461 : }
2462 0 : }
2463 :
2464 0 : SDGrid::CacheManager::~CacheManager()
2465 : {
2466 0 : if (onDuty) {
2467 0 : cache.close();
2468 : }
2469 0 : }
2470 :
2471 0 : SDGrid::CacheWriter::CacheWriter(Cache& cacheIn,
2472 0 : Bool onDutyIn)
2473 0 : : cache {cacheIn},
2474 0 : onDuty {onDutyIn}
2475 0 : {}
2476 :
2477 0 : SDGrid::CacheWriter::~CacheWriter()
2478 : {
2479 0 : if (onDuty) {
2480 0 : cache.storeRowPixel();
2481 : }
2482 0 : }
2483 :
2484 0 : const String & SDGrid::toString(const ConvertFirst convertFirst) {
2485 : static const std::array<String,3> name {
2486 : "NEVER",
2487 : "ALWAYS",
2488 : "AUTO"
2489 0 : };
2490 0 : switch(convertFirst){
2491 0 : case ConvertFirst::NEVER:
2492 : case ConvertFirst::ALWAYS:
2493 : case ConvertFirst::AUTO:
2494 0 : return name[static_cast<size_t>(convertFirst)];
2495 0 : default:
2496 0 : String errMsg {"Illegal ConvertFirst enum: "};
2497 0 : errMsg += String::toString(static_cast<Int>(convertFirst));
2498 0 : throw AipsError(
2499 : errMsg,
2500 : __FILE__,
2501 : __LINE__,
2502 : AipsError::Category::INVALID_ARGUMENT
2503 0 : );
2504 : // Avoid potential compiler warning
2505 : return name[static_cast<size_t>(ConvertFirst::NEVER)];
2506 : }
2507 : }
2508 :
2509 0 : SDGrid::ConvertFirst SDGrid::fromString(const String & name) {
2510 : static const std::array<ConvertFirst,3> schemes {
2511 : ConvertFirst::NEVER,
2512 : ConvertFirst::ALWAYS,
2513 : ConvertFirst::AUTO
2514 : };
2515 0 : for ( const auto scheme : schemes ) {
2516 0 : if (name == toString(scheme)) return scheme;
2517 : }
2518 0 : String errMsg {"Illegal ConvertFirst name: "};
2519 0 : errMsg += name;
2520 0 : throw AipsError(
2521 : errMsg,
2522 : __FILE__,
2523 : __LINE__,
2524 : AipsError::Category::INVALID_ARGUMENT
2525 0 : );
2526 : // Avoid potential compiler warning
2527 : return ConvertFirst::NEVER;
2528 0 : }
2529 :
2530 0 : void SDGrid::setConvertFirst(const casacore::String &name) {
2531 0 : processingScheme = fromString(name);
2532 0 : }
2533 :
2534 0 : Bool SDGrid::mustConvertPointingColumn(
2535 : const MeasurementSet &ms
2536 : ) {
2537 : const auto haveCachedSpectraPixelCoordinates =
2538 0 : cacheIsEnabled and cache.isReadable();
2539 0 : if (haveCachedSpectraPixelCoordinates) return False;
2540 :
2541 0 : switch(processingScheme){
2542 0 : case ConvertFirst::ALWAYS: return True;
2543 0 : case ConvertFirst::NEVER: return False;
2544 0 : case ConvertFirst::AUTO:
2545 : {
2546 0 : const auto nPointings = ms.pointing().nrow();
2547 0 : const auto nSelectedDataRows = ms.nrow();
2548 0 : return nSelectedDataRows > nPointings ? True : False;
2549 : }
2550 0 : default:
2551 0 : String errMsg {"Unexpected invalid state: "};
2552 0 : errMsg += "ConvertFirst processingScheme=";
2553 0 : errMsg += String::toString<Int>(static_cast<Int>(processingScheme));
2554 0 : errMsg += " ms=" + ms.tableName();
2555 0 : throw AipsError(
2556 : errMsg,
2557 : __FILE__,
2558 : __LINE__,
2559 : AipsError::Category::GENERAL
2560 0 : );
2561 : }
2562 : // Avoid potential compiler warning
2563 : return False;
2564 : }
2565 :
2566 0 : void SDGrid::handleNewMs(
2567 : ROVisibilityIterator &vi,
2568 : const ImageInterface<Complex>& image) {
2569 :
2570 : // Synchronize spatial coordinates cache
2571 0 : if (cacheIsEnabled) cache.newMS(vi.ms());
2572 :
2573 : // Handle interpolate-convert processing scheme
2574 0 : if (mustConvertPointingColumn(vi.ms())) {
2575 0 : const auto columnEnum = MSPointing::columnType(pointingDirCol_p);
2576 : const auto refType =
2577 0 : image.coordinates().directionCoordinate().directionType();
2578 0 : convertPointingColumn(vi.ms(), columnEnum, refType);
2579 : }
2580 : else {
2581 0 : ramPointingTable = MSPointing();
2582 0 : ramPointingColumnsPtr.reset();
2583 : }
2584 0 : }
2585 :
2586 : namespace convert_pointing_column_helpers {
2587 : using DirectionComputer = std::function<MDirection (rownr_t)>;
2588 :
2589 : DirectionComputer
2590 0 : metaDirectionComputer(
2591 : const MSPointingColumns &pointingColumns,
2592 : MSPointing::PredefinedColumns columnEnum) {
2593 : using std::placeholders::_1;
2594 : using Column = MSPointing::PredefinedColumns;
2595 0 : const auto & encoderDirections = pointingColumns.encoderMeas();
2596 0 : switch(columnEnum) {
2597 0 : case Column::DIRECTION:
2598 0 : return std::bind(
2599 0 : &MSPointingColumns::directionMeas,
2600 0 : &pointingColumns,
2601 : _1,
2602 0 : Double {0.0}
2603 0 : );
2604 : break;
2605 0 : case Column::TARGET:
2606 0 : return std::bind(
2607 0 : &MSPointingColumns::targetMeas,
2608 0 : &pointingColumns,
2609 : _1,
2610 0 : Double {0.0}
2611 0 : );
2612 : break;
2613 0 : case Column::SOURCE_OFFSET:
2614 0 : return std::bind(
2615 0 : &MSPointingColumns::sourceOffsetMeas,
2616 0 : &pointingColumns,
2617 : _1,
2618 0 : Double {0.0}
2619 0 : );
2620 : break;
2621 0 : case Column::POINTING_OFFSET:
2622 0 : return std::bind(
2623 0 : &MSPointingColumns::pointingOffsetMeas,
2624 0 : &pointingColumns,
2625 : _1,
2626 0 : Double {0.0}
2627 0 : );
2628 : break;
2629 0 : case Column::ENCODER:
2630 0 : return std::bind(
2631 0 : &ScalarMeasColumn<MDirection>::operator(),
2632 0 : &encoderDirections,
2633 : _1
2634 0 : );
2635 : break;
2636 0 : default:
2637 0 : throw AipsError(
2638 0 : String("Illegal Pointing Column Enum: " + String::toString(columnEnum)),
2639 : AipsError::INVALID_ARGUMENT
2640 0 : );
2641 : }
2642 : }
2643 :
2644 : // DirectionArchiver
2645 : class DirectionArchiver {
2646 : public:
2647 : virtual void put(rownr_t row, const MDirection & dir) = 0;
2648 : };
2649 :
2650 : // Derived Templated Class,
2651 : // implementing the shared constructor.
2652 : template<typename ColumnType, typename CellType>
2653 : class DirectionArchiver_ : public DirectionArchiver {
2654 : public:
2655 0 : DirectionArchiver_(ColumnType &columnIn)
2656 0 : : column {columnIn}
2657 0 : {}
2658 : void put(rownr_t row, const MDirection & dir);
2659 : private:
2660 : ColumnType & column;
2661 : };
2662 :
2663 : // Specializations of "put" member
2664 : template<>
2665 0 : void DirectionArchiver_<MDirection::ScalarColumn, MDirection>::put(
2666 : rownr_t row, const MDirection &dir) {
2667 0 : column.put(row, dir);
2668 0 : }
2669 :
2670 : template<>
2671 0 : void DirectionArchiver_<MDirection::ArrayColumn, Array<MDirection>>::put(
2672 : rownr_t row, const MDirection &dir) {
2673 0 : column.put(row, Vector<MDirection> {dir});
2674 0 : }
2675 :
2676 : // Template instantiations for the types we are interested in.
2677 : template class DirectionArchiver_ <MDirection::ArrayColumn, Array<MDirection>>;
2678 : template class DirectionArchiver_ <MDirection::ScalarColumn, MDirection>;
2679 :
2680 : // Aliases
2681 : using ScalarArchiver =
2682 : DirectionArchiver_ <MDirection::ScalarColumn, MDirection>;
2683 : using ArrayArchiver =
2684 : DirectionArchiver_ <MDirection::ArrayColumn, Array<MDirection>>;
2685 :
2686 :
2687 : struct ArchiverFactory {
2688 : static DirectionArchiver * createArchiver(
2689 : TableMeasColumn & column
2690 : );
2691 : };
2692 :
2693 : DirectionArchiver *
2694 0 : ArchiverFactory::createArchiver(TableMeasColumn &column) {
2695 : try {
2696 0 : auto & scalarColumn = dynamic_cast<MDirection::ScalarColumn &>(column);
2697 0 : return new ScalarArchiver(scalarColumn);
2698 : }
2699 0 : catch(std::bad_cast & exception) {
2700 0 : auto & arrayColumn = dynamic_cast<MDirection::ArrayColumn &>(column);
2701 0 : return new ArrayArchiver(arrayColumn);
2702 0 : }
2703 : }
2704 :
2705 : TableMeasColumn &
2706 0 : columnData(
2707 : MSPointingColumns & pointingColumns,
2708 : MSPointing::PredefinedColumns columnEnum) {
2709 : using Column = MSPointing::PredefinedColumns;
2710 0 : switch(columnEnum){
2711 : // Array Columns
2712 0 : case Column::DIRECTION:
2713 0 : return pointingColumns.directionMeasCol();
2714 0 : case Column::TARGET:
2715 0 : return pointingColumns.targetMeasCol();
2716 0 : case Column::SOURCE_OFFSET:
2717 0 : return pointingColumns.sourceOffsetMeasCol();
2718 0 : case Column::POINTING_OFFSET:
2719 0 : return pointingColumns.pointingOffsetMeasCol();
2720 : // Scalar Column
2721 0 : case Column::ENCODER:
2722 0 : return pointingColumns.encoderMeas();
2723 0 : default:
2724 : {
2725 0 : LogIO logger {LogOrigin {"columnData"} };
2726 : logger << LogIO::EXCEPTION
2727 : << "Expected a column of directions, got: " << MSPointing::columnName(columnEnum)
2728 0 : << LogIO::POST;
2729 :
2730 : // This is just to silence the following compiler warning:
2731 : // warning: control reaches end of non-void function [-Wreturn-type]
2732 0 : return pointingColumns.directionMeasCol();
2733 0 : }
2734 : }
2735 : }
2736 :
2737 : } // namespace convert_pointing_column_helpers
2738 : using namespace convert_pointing_column_helpers;
2739 :
2740 0 : void SDGrid::convertPointingColumn(
2741 : const MeasurementSet & ms,
2742 : const MSPointingEnums::PredefinedColumns columnEnum,
2743 : const MDirection::Types refType) {
2744 :
2745 0 : LogIO logger {LogOrigin {"SDGrid", "convertPointingColumn"}};
2746 0 : logger << "Start" << LogIO::POST;
2747 :
2748 0 : initRamPointingTable(ms.pointing(), columnEnum, refType);
2749 :
2750 : // Setup helper tools
2751 : // ---- Conversion Tools
2752 0 : MeasFrame pointingMeasurements;
2753 0 : MDirection::Convert convertToImageDirectionRef;
2754 0 : std::tie(
2755 : pointingMeasurements,
2756 : convertToImageDirectionRef
2757 0 : ) = setupConversionTools(ms, refType);
2758 :
2759 : // ---- Direction Computer
2760 0 : MSPointingColumns pointingColumns {ms.pointing()};
2761 : auto userSpecifiedDirection =
2762 0 : metaDirectionComputer(pointingColumns, columnEnum);
2763 :
2764 : // ---- Direction Archiver
2765 0 : MSPointingColumns ramPointingColumns {ramPointingTable};
2766 : std::unique_ptr<DirectionArchiver> correspondingRamPointingColumn {
2767 : ArchiverFactory::createArchiver(
2768 : columnData(
2769 : ramPointingColumns,
2770 : columnEnum
2771 : )
2772 : )
2773 0 : };
2774 :
2775 : // Convert directions stored in user-specified
2776 : // POINTING column (of some direction),
2777 : // and store the result into the corresponding column
2778 : // of the RAM POINTING table
2779 0 : const auto & epoch = pointingColumns.timeMeas();
2780 0 : const auto & antennaId = pointingColumns.antennaId();
2781 :
2782 0 : MSAntennaColumns antennaColumns(ms.antenna());
2783 0 : const auto & antennaPosition = antennaColumns.positionMeas();
2784 :
2785 : // Main loop control
2786 0 : const auto pointingRows = ms.pointing().nrow();
2787 0 : constexpr Int invalidAntennaId = -1;
2788 0 : auto previousPointing_AntennaId {invalidAntennaId};
2789 :
2790 : // Main loop
2791 0 : for (rownr_t pointingRow = 0; pointingRow < pointingRows ; ++pointingRow){
2792 : // Collect pointing measurements
2793 0 : const auto pointing_Epoch = epoch(pointingRow);
2794 0 : pointingMeasurements.resetEpoch(pointing_Epoch);
2795 :
2796 0 : const auto pointing_AntennaId = antennaId(pointingRow);
2797 0 : const auto antennaChanged =
2798 : ( pointing_AntennaId != previousPointing_AntennaId );
2799 0 : if (antennaChanged) {
2800 : const auto new_AntennaPosition =
2801 0 : antennaPosition(pointing_AntennaId);
2802 0 : pointingMeasurements.resetPosition(new_AntennaPosition);
2803 0 : previousPointing_AntennaId = pointing_AntennaId;
2804 0 : }
2805 :
2806 : // Convert
2807 : const auto convertedDirection =
2808 : convertToImageDirectionRef(
2809 0 : userSpecifiedDirection(pointingRow)
2810 0 : );
2811 :
2812 : // Store
2813 0 : correspondingRamPointingColumn->put(
2814 : pointingRow,
2815 : convertedDirection
2816 : );
2817 0 : }
2818 :
2819 0 : logger << "Done" << LogIO::POST;
2820 0 : }
2821 :
2822 0 : void SDGrid::initRamPointingTable(
2823 : const MSPointing & pointingTable,
2824 : const MSPointingEnums::PredefinedColumns columnEnum,
2825 : const MDirection::Types refType) {
2826 :
2827 0 : LogIO logger {LogOrigin {"SDGrid", "initRamPointingTable"}};
2828 0 : logger << "Start" << LogIO::POST;
2829 :
2830 0 : constexpr auto doNotCopyRows = True;
2831 0 : ramPointingTable = pointingTable.copyToMemoryTable(
2832 0 : pointingTable.tableName() +
2833 0 : "." + MSPointing::columnName(columnEnum) +
2834 0 : "." + MDirection::showType(refType),
2835 : doNotCopyRows
2836 0 : );
2837 0 : ramPointingColumnsPtr.reset( new MSPointingColumns {ramPointingTable});
2838 :
2839 0 : MSPointingColumns pointingColumns {pointingTable};
2840 :
2841 0 : ramPointingColumnsPtr->setDirectionRef(refType);
2842 0 : ramPointingColumnsPtr->setEncoderDirectionRef(refType);
2843 :
2844 0 : ramPointingTable.addRow(pointingTable.nrow());
2845 :
2846 0 : ramPointingColumnsPtr->antennaId().putColumn(pointingColumns.antennaId());
2847 0 : ramPointingColumnsPtr->time().putColumn(pointingColumns.time());
2848 0 : ramPointingColumnsPtr->interval().putColumn(pointingColumns.interval());
2849 0 : ramPointingColumnsPtr->numPoly().fillColumn(0);
2850 :
2851 0 : logger << "Done" << LogIO::POST;
2852 0 : }
2853 :
2854 : std::pair<MeasFrame,MDirection::Convert>
2855 0 : SDGrid::setupConversionTools(
2856 : const MeasurementSet & ms,
2857 : const casacore::MDirection::Types refType) {
2858 :
2859 0 : LogIO logger {LogOrigin {"SDGrid", "setupConversionTools"}};
2860 0 : logger << "Start" << LogIO::POST;
2861 :
2862 0 : MSPointingColumns pointingColumns(ms.pointing());
2863 0 : MSAntennaColumns antennaColumns(ms.antenna());
2864 :
2865 0 : auto firstPointing_Epoch = pointingColumns.timeMeas()(0);
2866 0 : auto firstPointing_AntennaId = pointingColumns.antennaId()(0);
2867 0 : auto firstPointing_AntennaPosition = antennaColumns.positionMeas()(
2868 : static_cast<rownr_t>(firstPointing_AntennaId)
2869 0 : );
2870 :
2871 0 : MeasFrame measFrame(firstPointing_Epoch, firstPointing_AntennaPosition);
2872 :
2873 : // Direction Conversion Machine
2874 0 : MDirection::Ref dstRef(refType, measFrame);
2875 0 : auto firstPointing_Direction = pointingColumns.directionMeas(0);
2876 0 : MDirection::Convert convert(firstPointing_Direction, dstRef);
2877 :
2878 : // ---- Perform 1 dummy conversion at a time far before
2879 : // the first pointing time, so that when we will next convert
2880 : // the first user-specified direction,
2881 : // we are sure that values cached in static variables
2882 : // of casacore functions like dUT1() will be cleared
2883 : static const MVEpoch oneYear {
2884 0 : Quantity {
2885 : 365,
2886 0 : Unit { "d" }
2887 : }
2888 0 : };
2889 : const MEpoch dummy_Epoch {
2890 0 : firstPointing_Epoch.getValue() - oneYear,
2891 0 : firstPointing_Epoch.getRef()
2892 0 : };
2893 0 : measFrame.resetEpoch(dummy_Epoch);
2894 0 : const auto dummy_Direction = convert();
2895 :
2896 0 : logger << "Done" << LogIO::POST;
2897 0 : return std::make_pair(measFrame, convert);
2898 0 : }
2899 :
2900 : } // casa namespace
|
