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