Line data Source code
1 : //# SDGrid.cc: Implementation of SDGrid class
2 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Arrays/Array.h>
29 : #include <casacore/casa/Arrays/ArrayLogical.h>
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Arrays/Cube.h>
32 : #include <casacore/casa/Arrays/MaskedArray.h>
33 : #include <casacore/casa/Arrays/Matrix.h>
34 : #include <casacore/casa/Arrays/MatrixIter.h>
35 : #include <casacore/casa/Arrays/Slice.h>
36 : #include <casacore/casa/Arrays/Vector.h>
37 : #include <casacore/casa/BasicSL/Constants.h>
38 : #include <casacore/casa/BasicSL/String.h>
39 : #include <casacore/casa/Containers/Block.h>
40 : #include <casacore/casa/Exceptions/Error.h>
41 : #include <casacore/casa/OS/Timer.h>
42 : #include <casacore/casa/Quanta/MVAngle.h>
43 : #include <casacore/casa/Quanta/MVTime.h>
44 : #include <casacore/casa/Quanta/UnitMap.h>
45 : #include <casacore/casa/Quanta/UnitVal.h>
46 : #include <sstream>
47 : #include <casacore/casa/Utilities/Assert.h>
48 :
49 : #include <components/ComponentModels/ConstantSpectrum.h>
50 : #include <components/ComponentModels/Flux.h>
51 : #include <components/ComponentModels/PointShape.h>
52 :
53 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
54 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
55 : #include <casacore/coordinates/Coordinates/Projection.h>
56 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
57 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
58 :
59 : #include <casacore/images/Images/ImageInterface.h>
60 : #include <casacore/images/Images/PagedImage.h>
61 : #include <casacore/images/Images/TempImage.h>
62 :
63 : #include <casacore/lattices/Lattices/ArrayLattice.h>
64 : #include <casacore/lattices/Lattices/LatticeCache.h>
65 : #include <casacore/lattices/Lattices/LatticeIterator.h>
66 : #include <casacore/lattices/Lattices/LatticeStepper.h>
67 : #include <casacore/lattices/Lattices/SubLattice.h>
68 : #include <casacore/lattices/LEL/LatticeExpr.h>
69 : #include <casacore/lattices/LRegions/LCBox.h>
70 :
71 : #include <casacore/measures/Measures/Stokes.h>
72 : #include <casacore/ms/MeasurementSets/MSColumns.h>
73 : #include <msvis/MSVis/StokesVector.h>
74 : #include <msvis/MSVis/VisBuffer2.h>
75 : #include <msvis/MSVis/VisibilityIterator2.h>
76 : #include <casacore/scimath/Mathematics/RigidVector.h>
77 : #include <synthesis/TransformMachines2/SDGrid.h>
78 : #include <synthesis/TransformMachines2/SkyJones.h>
79 : #include <synthesis/TransformMachines/StokesImageUtil.h>
80 :
81 : using namespace casacore;
82 : namespace casa {
83 : namespace refim {//# namespace for imaging refactor
84 :
85 0 : SDGrid::SDGrid(SkyJones& sj, Int icachesize, Int itilesize,
86 0 : String iconvType, Int userSupport, Bool useImagingWeight)
87 : : FTMachine(), sj_p(&sj), imageCache(0), wImageCache(0),
88 : cachesize(icachesize), tilesize(itilesize),
89 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
90 : pointingToImage(0), userSetSupport_p(userSupport),
91 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
92 : minWeight_p(0.), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
93 0 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(false)
94 : {
95 0 : lastIndex_p=0;
96 0 : }
97 :
98 36 : SDGrid::SDGrid(MPosition& mLocation, SkyJones& sj, Int icachesize, Int itilesize,
99 36 : String iconvType, Int userSupport, Float minweight, Bool clipminmax, Bool useImagingWeight)
100 : : FTMachine(), sj_p(&sj), imageCache(0), wImageCache(0),
101 : cachesize(icachesize), tilesize(itilesize),
102 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
103 : pointingToImage(0), userSetSupport_p(userSupport),
104 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
105 : minWeight_p(minweight), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
106 36 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(clipminmax)
107 : {
108 36 : mLocation_p=mLocation;
109 36 : lastIndex_p=0;
110 36 : }
111 :
112 0 : SDGrid::SDGrid(Int icachesize, Int itilesize,
113 0 : String iconvType, Int userSupport, Bool useImagingWeight)
114 : : FTMachine(), sj_p(0), imageCache(0), wImageCache(0),
115 : cachesize(icachesize), tilesize(itilesize),
116 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
117 : pointingToImage(0), userSetSupport_p(userSupport),
118 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
119 : minWeight_p(0.), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
120 0 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(false)
121 : {
122 0 : lastIndex_p=0;
123 0 : }
124 :
125 216 : SDGrid::SDGrid(MPosition &mLocation, Int icachesize, Int itilesize,
126 216 : String iconvType, Int userSupport, Float minweight, Bool clipminmax, Bool useImagingWeight)
127 : : FTMachine(), sj_p(0), imageCache(0), wImageCache(0),
128 : cachesize(icachesize), tilesize(itilesize),
129 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
130 : pointingToImage(0), userSetSupport_p(userSupport),
131 : truncate_p(-1.0), gwidth_p(0.0), jwidth_p(0.0),
132 : minWeight_p(minweight), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1),
133 : msId_p(-1),
134 216 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(clipminmax)
135 : {
136 216 : mLocation_p=mLocation;
137 216 : lastIndex_p=0;
138 216 : }
139 :
140 6 : SDGrid::SDGrid(MPosition &mLocation, Int icachesize, Int itilesize,
141 : String iconvType, Float truncate, Float gwidth, Float jwidth,
142 6 : Float minweight, Bool clipminmax, Bool useImagingWeight)
143 : : FTMachine(), sj_p(0), imageCache(0), wImageCache(0),
144 : cachesize(icachesize), tilesize(itilesize),
145 : isTiled(false), wImage(0), arrayLattice(0), wArrayLattice(0), lattice(0), wLattice(0), convType(iconvType),
146 : pointingToImage(0), userSetSupport_p(-1),
147 : truncate_p(truncate), gwidth_p(gwidth), jwidth_p(jwidth),
148 : minWeight_p(minweight), lastIndexPerAnt_p(), useImagingWeight_p(useImagingWeight), lastAntID_p(-1), msId_p(-1),
149 6 : isSplineInterpolationReady(false), interpolator(0), clipminmax_(clipminmax)
150 : {
151 6 : mLocation_p=mLocation;
152 6 : lastIndex_p=0;
153 6 : }
154 :
155 :
156 : //----------------------------------------------------------------------
157 0 : SDGrid& SDGrid::operator=(const SDGrid& other)
158 : {
159 0 : if(this!=&other) {
160 : //Do the base parameters
161 0 : FTMachine::operator=(other);
162 0 : sj_p=other.sj_p;
163 0 : imageCache=other.imageCache;
164 0 : wImage=other.wImage;
165 0 : wImageCache=other.wImageCache;
166 0 : cachesize=other.cachesize;
167 0 : tilesize=other.tilesize;
168 0 : isTiled=other.isTiled;
169 0 : lattice=other.lattice;
170 0 : arrayLattice=other.arrayLattice;
171 0 : wLattice=other.wLattice;
172 0 : wArrayLattice=other.wArrayLattice;
173 0 : convType=other.convType;
174 0 : if(other.wImage !=0)
175 0 : wImage=(TempImage<Float> *)other.wImage->cloneII();
176 : else
177 0 : wImage=0;
178 0 : pointingDirCol_p=other.pointingDirCol_p;
179 0 : pointingToImage=0;
180 0 : xyPos.resize();
181 0 : xyPos=other.xyPos;
182 0 : xyPosMovingOrig_p=other.xyPosMovingOrig_p;
183 0 : convFunc.resize();
184 0 : convFunc=other.convFunc;
185 0 : convSampling=other.convSampling;
186 0 : convSize=other.convSize;
187 0 : convSupport=other.convSupport;
188 0 : userSetSupport_p=other.userSetSupport_p;
189 0 : lastIndex_p=0;
190 0 : lastIndexPerAnt_p=0;
191 0 : lastAntID_p=-1;
192 0 : msId_p=-1;
193 0 : useImagingWeight_p=other.useImagingWeight_p;
194 0 : clipminmax_=other.clipminmax_;
195 : };
196 0 : return *this;
197 : };
198 :
199 516 : String SDGrid::name() const{
200 516 : return String("SDGrid");
201 : }
202 :
203 : //----------------------------------------------------------------------
204 : // Odds are that it changed.....
205 0 : Bool SDGrid::changed(const vi::VisBuffer2& /*vb*/) {
206 0 : return false;
207 : }
208 :
209 : //----------------------------------------------------------------------
210 0 : SDGrid::SDGrid(const SDGrid& other):FTMachine()
211 : {
212 0 : operator=(other);
213 0 : }
214 :
215 : #define NEED_UNDERSCORES
216 : #if defined(NEED_UNDERSCORES)
217 : #define grdsf grdsf_
218 : #define grdgauss grdgauss_
219 : #define grdjinc1 grdjinc1_
220 : #endif
221 :
222 : extern "C" {
223 : void grdsf(Double*, Double*);
224 : void grdgauss(Double*, Double*, Double*);
225 : void grdjinc1(Double*, Double*, Int*, Double*);
226 : }
227 :
228 : //----------------------------------------------------------------------
229 129 : void SDGrid::init() {
230 :
231 129 : logIO() << LogOrigin("SDGrid", "init") << LogIO::NORMAL;
232 :
233 : //pfile = fopen("ptdata.txt","w");
234 :
235 129 : ok();
236 :
237 : /*if((image->shape().product())>cachesize) {
238 : isTiled=true;
239 : }
240 : else {
241 : isTiled=false;
242 : }*/
243 129 : isTiled=false;
244 129 : nx = image->shape()(0);
245 129 : ny = image->shape()(1);
246 129 : npol = image->shape()(2);
247 129 : nchan = image->shape()(3);
248 :
249 129 : sumWeight.resize(npol, nchan);
250 :
251 : // Set up image cache needed for gridding. For BOX-car convolution
252 : // we can use non-overlapped tiles. Otherwise we need to use
253 : // overlapped tiles and additive gridding so that only increments
254 : // to a tile are written.
255 129 : if(imageCache) delete imageCache; imageCache=0;
256 :
257 129 : convType=downcase(convType);
258 129 : logIO() << "Convolution function : " << convType << LogIO::DEBUG1 << LogIO::POST;
259 129 : if(convType=="pb") {
260 : //cerr << "CNVFunc " << convFunc << endl;
261 :
262 : }
263 111 : else if(convType=="box") {
264 91 : convSupport=(userSetSupport_p >= 0) ? userSetSupport_p : 0;
265 91 : logIO() << "Support : " << convSupport << " pixels" << LogIO::POST;
266 91 : convSampling=100;
267 91 : convSize=convSampling*(2*convSupport+2);
268 91 : convFunc.resize(convSize);
269 91 : convFunc=0.0;
270 9191 : for (Int i=0;i<convSize/2;i++) {
271 9100 : convFunc(i)=1.0;
272 : }
273 : }
274 20 : else if(convType=="sf") {
275 : // SF
276 17 : convSupport=(userSetSupport_p >= 0) ? userSetSupport_p : 3;
277 17 : logIO() << "Support : " << convSupport << " pixels" << LogIO::POST;
278 17 : convSampling=100;
279 17 : convSize=convSampling*(2*convSupport+2);
280 17 : convFunc.resize(convSize);
281 17 : convFunc=0.0;
282 9917 : for (Int i=0;i<convSampling*convSupport;i++) {
283 9900 : Double nu=Double(i)/Double(convSupport*convSampling);
284 : Double val;
285 9900 : grdsf(&nu, &val);
286 9900 : convFunc(i)=(1.0-nu*nu)*val;
287 : }
288 : }
289 3 : else if(convType=="gauss") {
290 : // default is b=1.0 (Mangum et al. 2007)
291 1 : Double hwhm=(gwidth_p > 0.0) ? Double(gwidth_p) : sqrt(log(2.0));
292 1 : Float truncate=(truncate_p >= 0.0) ? truncate_p : 3.0*hwhm;
293 1 : convSampling=100;
294 1 : Int itruncate=(Int)(truncate*Double(convSampling)+0.5);
295 1 : logIO() << LogIO::DEBUG1 << "hwhm=" << hwhm << LogIO::POST;
296 1 : logIO() << LogIO::DEBUG1 << "itruncate=" << itruncate << LogIO::POST;
297 : //convSupport=(Int)(truncate+0.5);
298 1 : convSupport = (Int)(truncate);
299 1 : convSupport += (((truncate-(Float)convSupport) > 0.0) ? 1 : 0);
300 1 : convSize=convSampling*(2*convSupport+2);
301 1 : convFunc.resize(convSize);
302 1 : convFunc=0.0;
303 : Double val, x;
304 252 : for (Int i = 0 ; i <= itruncate ; i++) {
305 251 : x = Double(i)/Double(convSampling);
306 251 : grdgauss(&hwhm, &x, &val);
307 251 : convFunc(i)=val;
308 : }
309 :
310 : // String outfile = convType + ".dat";
311 : // ofstream ofs(outfile.c_str());
312 : // for (Int i = 0 ; i < convSize ; i++) {
313 : // ofs << i << " " << convFunc[i] << endl;
314 : // }
315 : // ofs.close();
316 : }
317 2 : else if (convType=="gjinc") {
318 : // default is b=2.52, c=1.55 (Mangum et al. 2007)
319 2 : Double hwhm = (gwidth_p > 0.0) ? Double(gwidth_p) : sqrt(log(2.0))*2.52;
320 2 : Double c = (jwidth_p > 0.0) ? Double(jwidth_p) : 1.55;
321 : //Float truncate = truncate_p;
322 2 : convSampling = 100;
323 2 : Int itruncate=(Int)(truncate_p*Double(convSampling)+0.5);
324 2 : logIO() << LogIO::DEBUG1 << "hwhm=" << hwhm << LogIO::POST;
325 2 : logIO() << LogIO::DEBUG1 << "c=" << c << LogIO::POST;
326 2 : logIO() << LogIO::DEBUG1 << "itruncate=" << itruncate << LogIO::POST;
327 : //convSupport=(truncate_p >= 0.0) ? (Int)(truncate_p+0.5) : (Int)(2*c+0.5);
328 2 : Float convSupportF = (truncate_p >= 0.0) ? truncate_p : (2*c);
329 2 : convSupport = (Int)convSupportF;
330 2 : convSupport += (((convSupportF-(Float)convSupport) > 0.0) ? 1 : 0);
331 2 : convSize=convSampling*(2*convSupport+2);
332 2 : convFunc.resize(convSize);
333 2 : convFunc=0.0;
334 : //UNUSED: Double r;
335 : Double x, val1, val2;
336 2 : Int normalize = 1;
337 2 : if (itruncate >= 0) {
338 0 : for (Int i = 0 ; i < itruncate ; i++) {
339 0 : x = Double(i) / Double(convSampling);
340 : //r = Double(i) / (Double(hwhm)*Double(convSampling));
341 0 : grdgauss(&hwhm, &x, &val1);
342 0 : grdjinc1(&c, &x, &normalize, &val2);
343 0 : convFunc(i) = val1 * val2;
344 : }
345 : }
346 : else {
347 : // default is to truncate at first null
348 382 : for (Int i=0;i<convSize;i++) {
349 382 : x = Double(i) / Double(convSampling);
350 : //r = Double(i) / (Double(hwhm)*Double(convSampling));
351 382 : grdjinc1(&c, &x, &normalize, &val2);
352 382 : if (val2 <= 0.0) {
353 2 : logIO() << LogIO::DEBUG1 << "convFunc is automatically truncated at radius " << x << LogIO::POST;
354 2 : break;
355 : }
356 380 : grdgauss(&hwhm, &x, &val1);
357 380 : convFunc(i) = val1 * val2;
358 : }
359 : }
360 :
361 : // String outfile = convType + ".dat";
362 : // ofstream ofs(outfile.c_str());
363 : // for (Int i = 0 ; i < convSize ; i++) {
364 : // ofs << i << " " << convFunc[i] << endl;
365 : // }
366 : // ofs.close();
367 : }
368 : else {
369 0 : logIO_p << "Unknown convolution function " << convType << LogIO::EXCEPTION;
370 : }
371 :
372 129 : if(wImage) delete wImage; wImage=0;
373 129 : wImage = new TempImage<Float>(image->shape(), image->coordinates());
374 :
375 : /*if(isTiled) {
376 : Float tileOverlap=0.5;
377 : if(convType=="box") {
378 : tileOverlap=0.0;
379 : }
380 : else {
381 : tileOverlap=0.5;
382 : tilesize=max(12,tilesize);
383 : }
384 : IPosition tileShape=IPosition(4,tilesize,tilesize,npol,nchan);
385 : Vector<Float> tileOverlapVec(4);
386 : tileOverlapVec=0.0;
387 : tileOverlapVec(0)=tileOverlap;
388 : tileOverlapVec(1)=tileOverlap;
389 : imageCache=new LatticeCache <Complex> (*image, cachesize, tileShape,
390 : tileOverlapVec,
391 : (tileOverlap>0.0));
392 :
393 : wImageCache=new LatticeCache <Float> (*wImage, cachesize, tileShape,
394 : tileOverlapVec,
395 : (tileOverlap>0.0));
396 :
397 : }
398 : */
399 129 : }
400 :
401 : // This is nasty, we should use CountedPointers here.
402 516 : SDGrid::~SDGrid() {
403 : //fclose(pfile);
404 258 : if (imageCache) delete imageCache; imageCache = 0;
405 258 : if (arrayLattice) delete arrayLattice; arrayLattice = 0;
406 258 : if (wImage) delete wImage; wImage = 0;
407 258 : if (wImageCache) delete wImageCache; wImageCache = 0;
408 258 : if (wArrayLattice) delete wArrayLattice; wArrayLattice = 0;
409 258 : if (interpolator) delete interpolator; interpolator = 0;
410 516 : }
411 :
412 18 : void SDGrid::findPBAsConvFunction(const ImageInterface<Complex>& image,
413 : const vi::VisBuffer2& vb) {
414 :
415 : // Get the coordinate system and increase the sampling by
416 : // a factor of ~ 100.
417 36 : CoordinateSystem coords(image.coordinates());
418 :
419 : // Set up the convolution function: make the buffer plenty
420 : // big so that we can trim it back
421 18 : convSupport=max(128, sj_p->support(vb, coords));
422 :
423 18 : convSampling=100;
424 18 : convSize=convSampling*convSupport;
425 :
426 : // Make a one dimensional image to calculate the
427 : // primary beam. We oversample this by a factor of
428 : // convSampling.
429 18 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
430 18 : AlwaysAssert(directionIndex>=0, AipsError);
431 36 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
432 36 : Vector<Double> sampling;
433 18 : sampling = dc.increment();
434 18 : sampling*=1.0/Double(convSampling);
435 18 : dc.setIncrement(sampling);
436 :
437 : // Set the reference value to the first pointing in the coordinate
438 : // system used in the POINTING table.
439 : {
440 18 : uInt row = 0;
441 :
442 : // reset lastAntID_p to use correct antenna position
443 18 : lastAntID_p = -1;
444 :
445 18 : const MSPointingColumns& act_mspc = vb.subtableColumns().pointing();
446 18 : Bool nullPointingTable = (act_mspc.nrow() < 1);
447 18 : Int pointIndex = -1;
448 18 : if (!nullPointingTable){
449 : //if(vb.newMS()) This thing is buggy...using msId change
450 18 : if (vb.msId() != msId_p) {
451 18 : lastIndex_p=0;
452 18 : if (lastIndexPerAnt_p.nelements() < (size_t)vb.nAntennas()) {
453 18 : lastIndexPerAnt_p.resize(vb.nAntennas());
454 : }
455 18 : lastIndexPerAnt_p = 0;
456 18 : msId_p = vb.msId();
457 : }
458 18 : pointIndex = getIndex(act_mspc, vb.time()(row), -1.0, vb.antenna1()(row));
459 18 : if (pointIndex < 0)
460 0 : pointIndex = getIndex(act_mspc, vb.time()(row), vb.timeInterval()(row), vb.antenna1()(row));
461 : }
462 18 : if (!nullPointingTable && ((pointIndex < 0) || (pointIndex >= Int(act_mspc.time().nrow())))) {
463 0 : ostringstream o;
464 0 : o << "Failed to find pointing information for time " <<
465 0 : MVTime(vb.time()(row)/86400.0);
466 0 : logIO_p << String(o) << LogIO::EXCEPTION;
467 : }
468 :
469 36 : MEpoch epoch(Quantity(vb.time()(row), "s"));
470 18 : if (!pointingToImage) {
471 18 : lastAntID_p = vb.antenna1()(row);
472 36 : MPosition pos = vb.subtableColumns().antenna().positionMeas()(lastAntID_p);
473 : //mFrame_p = MeasFrame(epoch, pos);
474 18 : (!mFrame_p.epoch()) ? mFrame_p.set(epoch) : mFrame_p.resetEpoch(epoch);
475 18 : (!mFrame_p.position()) ? mFrame_p.set(pos) : mFrame_p.resetPosition(pos);
476 18 : if (!nullPointingTable) {
477 18 : worldPosMeas = directionMeas(act_mspc, pointIndex);
478 : } else {
479 0 : worldPosMeas = vb.direction1()(row);
480 : }
481 :
482 : // Make a machine to convert from the worldPosMeas to the output
483 : // Direction Measure type for the relevant frame
484 36 : MDirection::Ref outRef(dc.directionType(), mFrame_p);
485 18 : pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
486 18 : if (!pointingToImage) {
487 0 : logIO_p << "Cannot make direction conversion machine" << LogIO::EXCEPTION;
488 : }
489 :
490 : } else {
491 0 : mFrame_p.resetEpoch(epoch);
492 0 : if (lastAntID_p != vb.antenna1()(row)) {
493 0 : logIO_p << LogIO::DEBUGGING
494 : << "updating antenna position. MS ID " << msId_p
495 : << ", last antenna ID " << lastAntID_p
496 0 : << " new antenna ID " << vb.antenna1()(row) << LogIO::POST;
497 0 : lastAntID_p = vb.antenna1()(row);
498 0 : MPosition pos = vb.subtableColumns().antenna().positionMeas()(lastAntID_p);
499 0 : mFrame_p.resetPosition(pos);
500 : }
501 : }
502 :
503 18 : if (!nullPointingTable) {
504 18 : worldPosMeas = (*pointingToImage)(directionMeas(act_mspc, pointIndex));
505 : } else {
506 0 : worldPosMeas = (*pointingToImage)(vb.direction1()(row));
507 : }
508 18 : delete pointingToImage;
509 18 : pointingToImage = 0;
510 : }
511 :
512 18 : directionCoord = coords.directionCoordinate(directionIndex);
513 : //make sure we use the same units
514 18 : worldPosMeas.set(dc.worldAxisUnits()(0));
515 :
516 : // Reference pixel may be modified in dc.setReferenceValue when
517 : // projection type is SFL. To take into account this effect,
518 : // keep original reference pixel here and subtract it from
519 : // the reference pixel after dc.setReferenceValue instead
520 : // of setting reference pixel to (0,0).
521 36 : Vector<Double> const originalReferencePixel = dc.referencePixel();
522 18 : dc.setReferenceValue(worldPosMeas.getAngle().getValue());
523 : //Vector<Double> unitVec(2);
524 : //unitVec=0.0;
525 : //dc.setReferencePixel(unitVec);
526 54 : Vector<Double> updatedReferencePixel = dc.referencePixel() - originalReferencePixel;
527 18 : dc.setReferencePixel(updatedReferencePixel);
528 :
529 18 : coords.replaceCoordinate(dc, directionIndex);
530 :
531 36 : IPosition pbShape(4, convSize, 2, 1, 1);
532 36 : IPosition start(4, 0, 0, 0, 0);
533 :
534 36 : TempImage<Complex> onedPB(pbShape, coords);
535 :
536 18 : onedPB.set(Complex(1.0, 0.0));
537 :
538 18 : AlwaysAssert(sj_p, AipsError);
539 18 : sj_p->apply(onedPB, onedPB, vb, 0);
540 :
541 36 : IPosition pbSlice(4, convSize, 1, 1, 1);
542 54 : Vector<Float> tempConvFunc=real(onedPB.getSlice(start, pbSlice, true));
543 : // Find number of significant points
544 18 : uInt cfLen=0;
545 7475 : for(uInt i=0;i<tempConvFunc.nelements();++i) {
546 7475 : if(tempConvFunc(i)<1e-3) break;
547 7457 : ++cfLen;
548 : }
549 18 : if(cfLen<1) {
550 : logIO() << LogIO::SEVERE
551 : << "Possible problem in primary beam calculation: no points in gridding function"
552 0 : << " - no points to be gridded on this image?" << LogIO::POST;
553 0 : cfLen=1;
554 : }
555 18 : Vector<Float> trimConvFunc=tempConvFunc(Slice(0,cfLen-1,1));
556 :
557 : // Now fill in the convolution function vector
558 18 : convSupport=cfLen/convSampling;
559 18 : convSize=convSampling*(2*convSupport+2);
560 18 : convFunc.resize(convSize);
561 18 : convFunc=0.0;
562 18 : convFunc(Slice(0,cfLen-1,1))=trimConvFunc(Slice(0,cfLen-1,1));
563 :
564 :
565 18 : }
566 :
567 : // Initialize for a transform from the Sky domain. This means that
568 : // we grid-correct, and FFT the image
569 0 : void SDGrid::initializeToVis(ImageInterface<Complex>& iimage,
570 : const vi::VisBuffer2& vb)
571 : {
572 0 : image=&iimage;
573 :
574 0 : ok();
575 :
576 0 : init();
577 :
578 0 : if(convType=="pb") {
579 0 : findPBAsConvFunction(*image, vb);
580 : }
581 :
582 : // reset msId_p and lastAntID_p to -1
583 : // this is to ensure correct antenna position in getXYPos
584 0 : msId_p = -1;
585 0 : lastAntID_p = -1;
586 :
587 : // Initialize the maps for polarization and channel. These maps
588 : // translate visibility indices into image indices
589 0 : initMaps(vb);
590 :
591 : // First get the CoordinateSystem for the image and then find
592 : // the DirectionCoordinate
593 0 : CoordinateSystem coords=image->coordinates();
594 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
595 0 : AlwaysAssert(directionIndex>=0, AipsError);
596 0 : directionCoord=coords.directionCoordinate(directionIndex);
597 : /*if((image->shape().product())>cachesize) {
598 : isTiled=true;
599 : }
600 : else {
601 : isTiled=false;
602 : }*/
603 0 : isTiled=false;
604 0 : nx = image->shape()(0);
605 0 : ny = image->shape()(1);
606 0 : npol = image->shape()(2);
607 0 : nchan = image->shape()(3);
608 :
609 : // If we are memory-based then read the image in and create an
610 : // ArrayLattice otherwise just use the PagedImage
611 : /*if(isTiled) {
612 : lattice=image;
613 : wLattice=wImage;
614 : }
615 : else*/
616 : {
617 : // Make the grid the correct shape and turn it into an array lattice
618 0 : IPosition gridShape(4, nx, ny, npol, nchan);
619 0 : griddedData.resize(gridShape);
620 0 : griddedData = Complex(0.0);
621 :
622 0 : wGriddedData.resize(gridShape);
623 0 : wGriddedData = 0.0;
624 :
625 0 : if(arrayLattice) delete arrayLattice; arrayLattice=0;
626 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
627 :
628 0 : if(wArrayLattice) delete wArrayLattice; wArrayLattice=0;
629 0 : wArrayLattice = new ArrayLattice<Float>(wGriddedData);
630 0 : wArrayLattice->set(0.0);
631 0 : wLattice=wArrayLattice;
632 :
633 : // Now find the SubLattice corresponding to the image
634 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2, (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
635 0 : IPosition stride(4, 1);
636 0 : IPosition trc(blc+image->shape()-stride);
637 0 : LCBox gridBox(blc, trc, gridShape);
638 0 : SubLattice<Complex> gridSub(*arrayLattice, gridBox, true);
639 :
640 : // Do the copy
641 0 : gridSub.copyData(*image);
642 :
643 0 : lattice=arrayLattice;
644 : }
645 0 : AlwaysAssert(lattice, AipsError);
646 0 : AlwaysAssert(wLattice, AipsError);
647 0 : }
648 :
649 0 : void SDGrid::finalizeToVis()
650 : {
651 : /*if(isTiled) {
652 :
653 : logIO() << LogOrigin("SDGrid", "finalizeToVis") << LogIO::NORMAL;
654 :
655 : AlwaysAssert(imageCache, AipsError);
656 : AlwaysAssert(image, AipsError);
657 : ostringstream o;
658 : imageCache->flush();
659 : imageCache->showCacheStatistics(o);
660 : logIO() << o.str() << LogIO::POST;
661 : }*/
662 0 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
663 0 : }
664 :
665 :
666 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
667 : // grid.
668 129 : void SDGrid::initializeToSky(ImageInterface<Complex>& iimage,
669 : Matrix<Float>& weight, const vi::VisBuffer2& vb)
670 : {
671 : // image always points to the image
672 129 : image=&iimage;
673 :
674 129 : ok();
675 :
676 129 : init();
677 :
678 129 : if(convType=="pb") {
679 18 : findPBAsConvFunction(*image, vb);
680 : }
681 :
682 : // reset msId_p and lastAntID_p to -1
683 : // this is to ensure correct antenna position in getXYPos
684 129 : msId_p = -1;
685 129 : lastAntID_p = -1;
686 :
687 : // Initialize the maps for polarization and channel. These maps
688 : // translate visibility indices into image indices
689 129 : initMaps(vb);
690 : //cerr << "ToSky cachesize " << cachesize << " im shape " << (image->shape().product()) << endl;
691 : /*if((image->shape().product())>cachesize) {
692 : isTiled=true;
693 : }
694 : else {
695 : isTiled=false;
696 : }
697 : */
698 : //////////////No longer using isTiled
699 129 : isTiled=false;
700 129 : nx = image->shape()(0);
701 129 : ny = image->shape()(1);
702 129 : npol = image->shape()(2);
703 129 : nchan = image->shape()(3);
704 :
705 129 : sumWeight=0.0;
706 129 : weight.resize(sumWeight.shape());
707 129 : weight=0.0;
708 :
709 : // First get the CoordinateSystem for the image and then find
710 : // the DirectionCoordinate
711 258 : CoordinateSystem coords=image->coordinates();
712 129 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
713 129 : AlwaysAssert(directionIndex>=0, AipsError);
714 129 : directionCoord=coords.directionCoordinate(directionIndex);
715 :
716 : // Initialize for in memory or to disk gridding. lattice will
717 : // point to the appropriate Lattice, either the ArrayLattice for
718 : // in memory gridding or to the image for to disk gridding.
719 : /*if(isTiled) {
720 : imageCache->flush();
721 : image->set(Complex(0.0));
722 : lattice=image;
723 : wLattice=wImage;
724 : }
725 : else*/
726 : {
727 258 : IPosition gridShape(4, nx, ny, npol, nchan);
728 258 : logIO() << LogOrigin("SDGrid", "initializeToSky", WHERE) << LogIO::DEBUGGING
729 258 : << "gridShape = " << gridShape << LogIO::POST;
730 129 : griddedData.resize(gridShape);
731 129 : griddedData=Complex(0.0);
732 129 : if(arrayLattice) delete arrayLattice; arrayLattice=0;
733 129 : arrayLattice = new ArrayLattice<Complex>(griddedData);
734 129 : lattice=arrayLattice;
735 129 : wGriddedData.resize(gridShape);
736 129 : wGriddedData=0.0;
737 129 : if(wArrayLattice) delete wArrayLattice; wArrayLattice=0;
738 129 : wArrayLattice = new ArrayLattice<Float>(wGriddedData);
739 129 : wLattice=wArrayLattice;
740 :
741 : // clipping related stuff
742 129 : if (clipminmax_) {
743 8 : gmin_.resize(gridShape);
744 8 : gmin_ = Complex(FLT_MAX);
745 8 : gmax_.resize(gridShape);
746 8 : gmax_ = Complex(-FLT_MAX);
747 8 : wmin_.resize(gridShape);
748 8 : wmin_ = 0.0f;
749 8 : wmax_.resize(gridShape);
750 8 : wmax_ = 0.0f;
751 8 : npoints_.resize(gridShape);
752 8 : npoints_ = 0;
753 : }
754 : }
755 129 : AlwaysAssert(lattice, AipsError);
756 129 : AlwaysAssert(wLattice, AipsError);
757 129 : }
758 :
759 129 : void SDGrid::finalizeToSky()
760 : {
761 :
762 : // Now we flush the cache and report statistics
763 : // For memory based, we don't write anything out yet.
764 : /*if(isTiled) {
765 : logIO() << LogOrigin("SDGrid", "finalizeToSky") << LogIO::NORMAL;
766 :
767 : AlwaysAssert(image, AipsError);
768 : AlwaysAssert(imageCache, AipsError);
769 : imageCache->flush();
770 : ostringstream o;
771 : imageCache->showCacheStatistics(o);
772 : logIO() << o.str() << LogIO::POST;
773 : }
774 : */
775 :
776 129 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
777 129 : }
778 :
779 0 : Array<Complex>* SDGrid::getDataPointer(const IPosition& centerLoc2D,
780 : Bool readonly) {
781 : Array<Complex>* result;
782 : // Is tiled: get tiles and set up offsets
783 0 : centerLoc(0)=centerLoc2D(0);
784 0 : centerLoc(1)=centerLoc2D(1);
785 0 : result=&imageCache->tile(offsetLoc, centerLoc, readonly);
786 0 : return result;
787 : }
788 0 : Array<Float>* SDGrid::getWDataPointer(const IPosition& centerLoc2D,
789 : Bool readonly) {
790 : Array<Float>* result;
791 : // Is tiled: get tiles and set up offsets
792 0 : centerLoc(0)=centerLoc2D(0);
793 0 : centerLoc(1)=centerLoc2D(1);
794 0 : result=&wImageCache->tile(offsetLoc, centerLoc, readonly);
795 0 : return result;
796 : }
797 :
798 : #define NEED_UNDERSCORES
799 : #if defined(NEED_UNDERSCORES)
800 : #define ggridsd ggridsd_
801 : #define dgridsd dgridsd_
802 : #define ggridsdclip ggridsdclip_
803 : #endif
804 :
805 : extern "C" {
806 : void ggridsd(Double*,
807 : const Complex*,
808 : Int*,
809 : Int*,
810 : Int*,
811 : const Int*,
812 : const Int*,
813 : const Float*,
814 : Int*,
815 : Int*,
816 : Complex*,
817 : Float*,
818 : Int*,
819 : Int*,
820 : Int *,
821 : Int *,
822 : Int*,
823 : Int*,
824 : Float*,
825 : Int*,
826 : Int*,
827 : Double*);
828 : void ggridsdclip(Double*,
829 : const Complex*,
830 : Int*,
831 : Int*,
832 : Int*,
833 : const Int*,
834 : const Int*,
835 : const Float*,
836 : Int*,
837 : Int*,
838 : Complex*,
839 : Float*,
840 : Int*,
841 : Complex*,
842 : Float*,
843 : Complex*,
844 : Float*,
845 : Int*,
846 : Int*,
847 : Int *,
848 : Int *,
849 : Int*,
850 : Int*,
851 : Float*,
852 : Int*,
853 : Int*,
854 : Double*);
855 : void dgridsd(Double*,
856 : Complex*,
857 : Int*,
858 : Int*,
859 : const Int*,
860 : const Int*,
861 : Int*,
862 : Int*,
863 : const Complex*,
864 : Int*,
865 : Int*,
866 : Int *,
867 : Int *,
868 : Int*,
869 : Int*,
870 : Float*,
871 : Int*,
872 : Int*);
873 : }
874 :
875 188425 : void SDGrid::put(const vi::VisBuffer2& vb, Int row, Bool dopsf,
876 : FTMachine::Type type)
877 : {
878 376850 : LogIO os(LogOrigin("SDGrid", "put"));
879 :
880 188425 : gridOk(convSupport);
881 :
882 : // There is no channel mapping cache in VI/VB2 version of FTMachine
883 : // Perform matchChannel everytime
884 188425 : matchChannel(vb);
885 :
886 : //No point in reading data if its not matching in frequency
887 188425 : if(max(chanMap)==-1)
888 19215 : return;
889 :
890 338420 : Matrix<Float> imagingweight;
891 : //imagingweight=&(vb.imagingWeight());
892 169210 : pickWeights(vb, imagingweight);
893 :
894 169210 : if(type==FTMachine::PSF || type==FTMachine::COVERAGE)
895 0 : dopsf=true;
896 169210 : if(dopsf) type=FTMachine::PSF;
897 338420 : Cube<Complex> data;
898 : //Fortran gridder need the flag as ints
899 338420 : Cube<Int> flags;
900 338420 : Matrix<Float> elWeight;
901 169210 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
902 : //cerr << "number of rows " << vb.nRow() << " data shape " << data.shape() << endl;
903 : Bool iswgtCopy;
904 : const Float *wgtStorage;
905 169210 : wgtStorage=elWeight.getStorage(iswgtCopy);
906 : Bool isCopy;
907 169210 : const Complex *datStorage=0;
908 169210 : if(!dopsf)
909 169210 : datStorage=data.getStorage(isCopy);
910 :
911 : // If row is -1 then we pass through all rows
912 : Int startRow, endRow, nRow;
913 169210 : if (row==-1) {
914 169210 : nRow=vb.nRows();
915 169210 : startRow=0;
916 169210 : endRow=nRow-1;
917 : } else {
918 0 : nRow=1;
919 0 : startRow=row;
920 0 : endRow=row;
921 : }
922 :
923 :
924 338420 : Vector<Int> rowFlags(vb.flagRow().nelements());
925 169210 : rowFlags=0;
926 345705 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
927 176495 : if(vb.flagRow()(rownr)) rowFlags(rownr)=1;
928 : }
929 :
930 : // Take care of translation of Bools to Integer
931 169210 : Int idopsf=0;
932 169210 : if(dopsf) idopsf=1;
933 :
934 : /*if(isTiled) {
935 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
936 :
937 : if(getXYPos(vb, rownr)) {
938 :
939 : IPosition centerLoc2D(2, Int(xyPos(0)), Int(xyPos(1)));
940 : Array<Complex>* dataPtr=getDataPointer(centerLoc2D, false);
941 : Array<Float>* wDataPtr=getWDataPointer(centerLoc2D, false);
942 : Int aNx=dataPtr->shape()(0);
943 : Int aNy=dataPtr->shape()(1);
944 : Vector<Double> actualPos(2);
945 : for (Int i=0;i<2;i++) {
946 : actualPos(i)=xyPos(i)-Double(offsetLoc(i));
947 : }
948 : // Now use FORTRAN to do the gridding. Remember to
949 : // ensure that the shape and offsets of the tile are
950 : // accounted for.
951 : {
952 : Bool del;
953 : // IPosition s(data.shape());
954 : const IPosition& fs=flags.shape();
955 : std::vector<Int> s(fs.begin(), fs.end());
956 :
957 : ggridsd(actualPos.getStorage(del),
958 : datStorage,
959 : &s[0],
960 : &s[1],
961 : &idopsf,
962 : flags.getStorage(del),
963 : rowFlags.getStorage(del),
964 : wgtStorage,
965 : &s[2],
966 : &rownr,
967 : dataPtr->getStorage(del),
968 : wDataPtr->getStorage(del),
969 : &aNx,
970 : &aNy,
971 : &npol,
972 : &nchan,
973 : &convSupport,
974 : &convSampling,
975 : convFunc.getStorage(del),
976 : chanMap.getStorage(del),
977 : polMap.getStorage(del),
978 : sumWeight.getStorage(del));
979 : }
980 : }
981 : }
982 : }
983 : else*/
984 : {
985 338420 : Matrix<Double> xyPositions(2, endRow-startRow+1);
986 169210 : xyPositions=-1e9; // make sure failed getXYPos does not fall on grid
987 345705 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
988 176495 : if(getXYPos(vb, rownr)) {
989 176495 : xyPositions(0, rownr)=xyPos(0);
990 176495 : xyPositions(1, rownr)=xyPos(1);
991 : }
992 : }
993 : {
994 : Bool del;
995 : // IPosition s(data.shape());
996 169210 : const IPosition& fs=flags.shape();
997 338420 : std::vector<Int> s(fs.begin(), fs.end());
998 : Bool datCopy, wgtCopy;
999 169210 : Complex * datStor=griddedData.getStorage(datCopy);
1000 169210 : Float * wgtStor=wGriddedData.getStorage(wgtCopy);
1001 :
1002 : //Bool call_ggridsd = !clipminmax_ || dopsf;
1003 169210 : Bool call_ggridsd = !clipminmax_;
1004 :
1005 169210 : if (call_ggridsd) {
1006 :
1007 338370 : ggridsd(xyPositions.getStorage(del),
1008 : datStorage,
1009 169185 : &s[0],
1010 169185 : &s[1],
1011 : &idopsf,
1012 169185 : flags.getStorage(del),
1013 169185 : rowFlags.getStorage(del),
1014 : wgtStorage,
1015 169185 : &s[2],
1016 : &row,
1017 : datStor,
1018 : wgtStor,
1019 : &nx,
1020 : &ny,
1021 : &npol,
1022 : &nchan,
1023 : &convSupport,
1024 : &convSampling,
1025 : convFunc.getStorage(del),
1026 : chanMap.getStorage(del),
1027 : polMap.getStorage(del),
1028 : sumWeight.getStorage(del));
1029 :
1030 : } else {
1031 : Bool gminCopy;
1032 25 : Complex *gminStor = gmin_.getStorage(gminCopy);
1033 : Bool gmaxCopy;
1034 25 : Complex *gmaxStor = gmax_.getStorage(gmaxCopy);
1035 : Bool wminCopy;
1036 25 : Float *wminStor = wmin_.getStorage(wminCopy);
1037 : Bool wmaxCopy;
1038 25 : Float *wmaxStor = wmax_.getStorage(wmaxCopy);
1039 : Bool npCopy;
1040 25 : Int *npStor = npoints_.getStorage(npCopy);
1041 :
1042 50 : ggridsdclip(xyPositions.getStorage(del),
1043 : datStorage,
1044 25 : &s[0],
1045 25 : &s[1],
1046 : &idopsf,
1047 25 : flags.getStorage(del),
1048 25 : rowFlags.getStorage(del),
1049 : wgtStorage,
1050 25 : &s[2],
1051 : &row,
1052 : datStor,
1053 : wgtStor,
1054 : npStor,
1055 : gminStor,
1056 : wminStor,
1057 : gmaxStor,
1058 : wmaxStor,
1059 : &nx,
1060 : &ny,
1061 : &npol,
1062 : &nchan,
1063 : &convSupport,
1064 : &convSampling,
1065 : convFunc.getStorage(del),
1066 : chanMap.getStorage(del),
1067 : polMap.getStorage(del),
1068 : sumWeight.getStorage(del));
1069 :
1070 25 : gmin_.putStorage(gminStor, gminCopy);
1071 25 : gmax_.putStorage(gmaxStor, gmaxCopy);
1072 25 : wmin_.putStorage(wminStor, wminCopy);
1073 25 : wmax_.putStorage(wmaxStor, wmaxCopy);
1074 25 : npoints_.putStorage(npStor, npCopy);
1075 : }
1076 169210 : griddedData.putStorage(datStor, datCopy);
1077 169210 : wGriddedData.putStorage(wgtStor, wgtCopy);
1078 : }
1079 : }
1080 169210 : if(!dopsf)
1081 169210 : data.freeStorage(datStorage, isCopy);
1082 169210 : elWeight.freeStorage(wgtStorage,iswgtCopy);
1083 :
1084 : }
1085 :
1086 0 : void SDGrid::get(vi::VisBuffer2& vb, Int row)
1087 : {
1088 0 : LogIO os(LogOrigin("SDGrid", "get"));
1089 :
1090 0 : gridOk(convSupport);
1091 :
1092 : // If row is -1 then we pass through all rows
1093 : Int startRow, endRow, nRow;
1094 0 : if (row==-1) {
1095 0 : nRow=vb.nRows();
1096 0 : startRow=0;
1097 0 : endRow=nRow-1;
1098 : // TODO: ask imager guru if commenting out the following line
1099 : // is safe for SDGrid
1100 : //unnecessary zeroing
1101 : //vb.modelVisCube()=Complex(0.0,0.0);
1102 : } else {
1103 0 : nRow=1;
1104 0 : startRow=row;
1105 0 : endRow=row;
1106 : // TODO: ask imager guru if commenting out the following line
1107 : // is safe for SDGrid
1108 : //unnecessary zeroing
1109 : //vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
1110 : }
1111 :
1112 : // There is no channel mapping cache in VI/VB2 version of FTMachine
1113 : // Perform matchChannel everytime
1114 0 : matchChannel(vb);
1115 :
1116 : //No point in reading data if its not matching in frequency
1117 0 : if(max(chanMap)==-1)
1118 0 : return;
1119 0 : Cube<Complex> data;
1120 0 : Cube<Int> flags;
1121 0 : getInterpolateArrays(vb, data, flags);
1122 :
1123 : Complex *datStorage;
1124 : Bool isCopy;
1125 0 : datStorage=data.getStorage(isCopy);
1126 : // NOTE: with MS V2.0 the pointing could change per antenna and timeslot
1127 : //
1128 0 : Vector<Int> rowFlags(vb.flagRow().nelements());
1129 0 : rowFlags=0;
1130 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1131 0 : if(vb.flagRow()(rownr)) rowFlags(rownr)=1;
1132 : //single dish yes ?
1133 0 : if(max(vb.uvw().column(rownr)) > 0.0) rowFlags(rownr)=1;
1134 : }
1135 :
1136 :
1137 : /*if(isTiled) {
1138 :
1139 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1140 :
1141 : if(getXYPos(vb, rownr)) {
1142 :
1143 : // Get the tile
1144 : IPosition centerLoc2D(2, Int(xyPos(0)), Int(xyPos(1)));
1145 : Array<Complex>* dataPtr=getDataPointer(centerLoc2D, true);
1146 : Int aNx=dataPtr->shape()(0);
1147 : Int aNy=dataPtr->shape()(1);
1148 :
1149 : // Now use FORTRAN to do the gridding. Remember to
1150 : // ensure that the shape and offsets of the tile are
1151 : // accounted for.
1152 : Bool del;
1153 : Vector<Double> actualPos(2);
1154 : for (Int i=0;i<2;i++) {
1155 : actualPos(i)=xyPos(i)-Double(offsetLoc(i));
1156 : }
1157 : // IPosition s(data.shape());
1158 : const IPosition& fs=data.shape();
1159 : std::vector<Int> s(fs.begin(), fs.end());
1160 :
1161 : dgridsd(actualPos.getStorage(del),
1162 : datStorage,
1163 : &s[0],
1164 : &s[1],
1165 : flags.getStorage(del),
1166 : rowFlags.getStorage(del),
1167 : &s[2],
1168 : &rownr,
1169 : dataPtr->getStorage(del),
1170 : &aNx,
1171 : &aNy,
1172 : &npol,
1173 : &nchan,
1174 : &convSupport,
1175 : &convSampling,
1176 : convFunc.getStorage(del),
1177 : chanMap.getStorage(del),
1178 : polMap.getStorage(del));
1179 : }
1180 : }
1181 : }
1182 : else*/
1183 : {
1184 0 : Matrix<Double> xyPositions(2, endRow-startRow+1);
1185 0 : xyPositions=-1e9;
1186 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
1187 0 : if(getXYPos(vb, rownr)) {
1188 0 : xyPositions(0, rownr)=xyPos(0);
1189 0 : xyPositions(1, rownr)=xyPos(1);
1190 : }
1191 : }
1192 :
1193 : Bool del;
1194 : // IPosition s(data.shape());
1195 0 : const IPosition& fs=data.shape();
1196 0 : std::vector<Int> s(fs.begin(), fs.end());
1197 0 : dgridsd(xyPositions.getStorage(del),
1198 : datStorage,
1199 0 : &s[0],
1200 0 : &s[1],
1201 0 : flags.getStorage(del),
1202 0 : rowFlags.getStorage(del),
1203 0 : &s[2],
1204 : &row,
1205 0 : griddedData.getStorage(del),
1206 : &nx,
1207 : &ny,
1208 : &npol,
1209 : &nchan,
1210 : &convSupport,
1211 : &convSampling,
1212 : convFunc.getStorage(del),
1213 : chanMap.getStorage(del),
1214 : polMap.getStorage(del));
1215 :
1216 0 : data.putStorage(datStorage, isCopy);
1217 : }
1218 0 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1219 : }
1220 :
1221 :
1222 :
1223 : // Make a plain straightforward honest-to-FSM image. This returns
1224 : // a complex image, without conversion to Stokes. The representation
1225 : // is that required for the visibilities.
1226 : //----------------------------------------------------------------------
1227 0 : void SDGrid::makeImage(FTMachine::Type type,
1228 : vi::VisibilityIterator2& vi,
1229 : ImageInterface<Complex>& theImage,
1230 : Matrix<Float>& weight) {
1231 :
1232 :
1233 0 : logIO() << LogOrigin("FTMachine", "makeImage0") << LogIO::NORMAL;
1234 :
1235 : // Loop over all visibilities and pixels
1236 0 : vi::VisBuffer2 *vb = vi.getVisBuffer();
1237 :
1238 : // Initialize put (i.e. transform to Sky) for this model
1239 0 : vi.origin();
1240 :
1241 0 : if(vb->polarizationFrame()==MSIter::Linear) {
1242 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
1243 : }
1244 : else {
1245 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
1246 : }
1247 0 : Bool useCorrected= !(MSMainColumns(vi.ms()).correctedData().isNull());
1248 0 : if((type==FTMachine::CORRECTED) && (!useCorrected))
1249 0 : type=FTMachine::OBSERVED;
1250 0 : Bool normalize=true;
1251 0 : if(type==FTMachine::COVERAGE)
1252 0 : normalize=false;
1253 :
1254 0 : Int Nx=theImage.shape()(0);
1255 0 : Int Ny=theImage.shape()(1);
1256 0 : Int Npol=theImage.shape()(2);
1257 0 : Int Nchan=theImage.shape()(3);
1258 0 : Double memtot=Double(HostInfo::memoryTotal(true))*1024.0; // return in kB
1259 0 : Int nchanInMem=Int(memtot/2.0/8.0/Double(Nx*Ny*Npol));
1260 0 : Int nloop=nchanInMem >= Nchan ? 1 : Nchan/nchanInMem+1;
1261 0 : ImageInterface<Complex> *imCopy=NULL;
1262 0 : IPosition blc(theImage.shape());
1263 0 : IPosition trc(theImage.shape());
1264 0 : blc-=blc; //set all values to 0
1265 0 : trc=theImage.shape();
1266 0 : trc-=1; // set trc to image size -1
1267 0 : if(nloop==1) {
1268 0 : imCopy=& theImage;
1269 0 : nchanInMem=Nchan;
1270 : }
1271 : else
1272 : logIO() << "Not enough memory to image in one go \n will process the image in "
1273 : << nloop
1274 : << " sections "
1275 0 : << LogIO::POST;
1276 :
1277 0 : weight.resize(Npol, Nchan);
1278 0 : Matrix<Float> wgtcopy(Npol, Nchan);
1279 :
1280 0 : Bool isWgtZero=true;
1281 0 : for (Int k=0; k < nloop; ++k){
1282 0 : Int bchan=k*nchanInMem;
1283 0 : Int echan=(k+1)*nchanInMem < Nchan ? (k+1)*nchanInMem-1 : Nchan-1;
1284 :
1285 0 : if(nloop > 1) {
1286 0 : blc[3]=bchan;
1287 0 : trc[3]=echan;
1288 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1289 0 : imCopy=new SubImage<Complex>(theImage, sl, true);
1290 0 : wgtcopy.resize(npol, echan-bchan+1);
1291 : }
1292 0 : vi.originChunks();
1293 0 : vi.origin();
1294 0 : initializeToSky(*imCopy,wgtcopy,*vb);
1295 :
1296 :
1297 : // for minmax clipping
1298 0 : logIO() << LogOrigin("SDGrid", "makeImage", WHERE) << LogIO::DEBUGGING
1299 0 : << "doclip_ = " << (clipminmax_ ? "TRUE" : "FALSE") << " (" << clipminmax_ << ")" << LogIO::POST;
1300 0 : if (clipminmax_) {
1301 0 : logIO() << LogOrigin("SDGRID", "makeImage", WHERE)
1302 0 : << LogIO::DEBUGGING << "use ggridsd2 for imaging" << LogIO::POST;
1303 : }
1304 :
1305 : // Loop over the visibilities, putting VisBuffers
1306 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1307 0 : for (vi.origin(); vi.more(); vi.next()) {
1308 :
1309 0 : switch(type) {
1310 0 : case FTMachine::RESIDUAL:
1311 0 : vb->setVisCube(vb->visCubeCorrected() - vb->visCubeModel());
1312 0 : put(*vb, -1, false);
1313 0 : break;
1314 0 : case FTMachine::MODEL:
1315 0 : put(*vb, -1, false, FTMachine::MODEL);
1316 0 : break;
1317 0 : case FTMachine::CORRECTED:
1318 0 : put(*vb, -1, false, FTMachine::CORRECTED);
1319 0 : break;
1320 0 : case FTMachine::PSF:
1321 0 : vb->setVisCube(Complex(1.0,0.0));
1322 0 : put(*vb, -1, true, FTMachine::PSF);
1323 0 : break;
1324 0 : case FTMachine::COVERAGE:
1325 0 : vb->setVisCube(Complex(1.0));
1326 0 : put(*vb, -1, true, FTMachine::COVERAGE);
1327 0 : break;
1328 0 : case FTMachine::OBSERVED:
1329 : default:
1330 0 : put(*vb, -1, false, FTMachine::OBSERVED);
1331 0 : break;
1332 : }
1333 : }
1334 : }
1335 0 : finalizeToSky();
1336 : // Normalize by dividing out weights, etc.
1337 0 : getImage(wgtcopy, normalize);
1338 0 : if(max(wgtcopy)==0.0){
1339 0 : if(nloop > 1)
1340 : logIO() << LogIO::WARN
1341 : << "No useful data in SDGrid: weights all zero for image slice " << k
1342 0 : << LogIO::POST;
1343 : }
1344 : else
1345 0 : isWgtZero=false;
1346 :
1347 0 : weight(Slice(0, Npol), Slice(bchan, echan-bchan+1))=wgtcopy;
1348 0 : if(nloop >1) delete imCopy;
1349 : }//loop k
1350 0 : if(isWgtZero)
1351 : logIO() << LogIO::SEVERE
1352 : << "No useful data in SDGrid: weights all zero"
1353 0 : << LogIO::POST;
1354 0 : }
1355 :
1356 :
1357 :
1358 : // Finalize : optionally normalize by weight image
1359 129 : ImageInterface<Complex>& SDGrid::getImage(Matrix<Float>& weights,
1360 : Bool normalize)
1361 : {
1362 129 : AlwaysAssert(lattice, AipsError);
1363 129 : AlwaysAssert(image, AipsError);
1364 :
1365 129 : logIO() << LogOrigin("SDGrid", "getImage") << LogIO::NORMAL;
1366 :
1367 : // execute minmax clipping
1368 129 : clipMinMax();
1369 :
1370 129 : weights.resize(sumWeight.shape());
1371 :
1372 129 : convertArray(weights,sumWeight);
1373 :
1374 : // If the weights are all zero then we cannot normalize
1375 : // otherwise we don't care.
1376 : ///////////////////////
1377 : /*{
1378 : PagedImage<Float> thisScreen(lattice->shape(), image->coordinates(), "TheData");
1379 : LatticeExpr<Float> le(abs(*lattice));
1380 : thisScreen.copyData(le);
1381 : PagedImage<Float> thisScreen2(lattice->shape(), image->coordinates(), "TheWeight");
1382 : LatticeExpr<Float> le2(abs(*wLattice));
1383 : thisScreen2.copyData(le2);
1384 : }*/
1385 : /////////////////////
1386 :
1387 129 : if(normalize) {
1388 0 : if(max(weights)==0.0) {
1389 : //logIO() << LogIO::WARN << "No useful data in SDGrid: weights all zero"
1390 : // << LogIO::POST;
1391 : //image->set(Complex(0.0));
1392 0 : return *image;
1393 : }
1394 : //Timer tim;
1395 : //tim.mark();
1396 : //lattice->copyData((LatticeExpr<Complex>)(iif((*wLattice<=minWeight_p), 0.0,
1397 : // (*lattice)/(*wLattice))));
1398 : //As we are not using disk based lattices anymore the above uses too much memory for
1399 : // ArrayLattice...it does not do a real inplace math but rather mutiple copies
1400 : // seem to be involved thus the less elegant but faster and less memory hog loop
1401 : // below.
1402 :
1403 0 : IPosition pos(4);
1404 0 : for (Int chan=0; chan < nchan; ++chan){
1405 0 : pos[3]=chan;
1406 0 : for( Int pol=0; pol < npol; ++ pol){
1407 0 : pos[2]=pol;
1408 0 : for (Int y=0; y < ny ; ++y){
1409 0 : pos[1]=y;
1410 0 : for (Int x=0; x < nx; ++x){
1411 0 : pos[0]=x;
1412 0 : Float wgt=wGriddedData(pos);
1413 0 : if(wgt > minWeight_p)
1414 0 : griddedData(pos)=griddedData(pos)/wgt;
1415 : else
1416 0 : griddedData(pos)=0.0;
1417 : }
1418 : }
1419 : }
1420 : }
1421 : //tim.show("After normalizing");
1422 : }
1423 :
1424 : //if(!isTiled)
1425 : {
1426 : // Now find the SubLattice corresponding to the image
1427 258 : IPosition gridShape(4, nx, ny, npol, nchan);
1428 258 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
1429 387 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1430 258 : IPosition stride(4, 1);
1431 387 : IPosition trc(blc+image->shape()-stride);
1432 258 : LCBox gridBox(blc, trc, gridShape);
1433 387 : SubLattice<Complex> gridSub(*arrayLattice, gridBox);
1434 :
1435 : // Do the copy
1436 129 : image->copyData(gridSub);
1437 : }
1438 :
1439 : // IMAGER MIGRATION
1440 : // set sumWeight to 1.0
1441 129 : weights = 1.0;
1442 :
1443 129 : return *image;
1444 : }
1445 :
1446 : // Return weights image
1447 129 : void SDGrid::getWeightImage(ImageInterface<Float>& weightImage, Matrix<Float>& weights)
1448 : {
1449 129 : AlwaysAssert(lattice, AipsError);
1450 129 : AlwaysAssert(image, AipsError);
1451 :
1452 129 : logIO() << LogOrigin("SDGrid", "getWeightImage") << LogIO::NORMAL;
1453 :
1454 129 : weights.resize(sumWeight.shape());
1455 : // IMAGER MIGRATION
1456 : // set sumWeight to 1.0
1457 129 : weights = 1.0;
1458 : // convertArray(weights,sumWeight);
1459 :
1460 129 : weightImage.copyData(*wArrayLattice);
1461 129 : }
1462 :
1463 258 : void SDGrid::ok() {
1464 258 : AlwaysAssert(image, AipsError);
1465 258 : }
1466 :
1467 : // Get the index into the pointing table for this time. Note that the
1468 : // in the pointing table, TIME specifies the beginning of the spanned
1469 : // time range, whereas for the main table, TIME is the centroid.
1470 : // Note that the behavior for multiple matches is not specified! i.e.
1471 : // if there are multiple matches, the index returned depends on the
1472 : // history of previous matches. It is deterministic but not obvious.
1473 : // One could cure this by searching but it would be considerably
1474 : // costlier.
1475 176513 : Int SDGrid::getIndex(const MSPointingColumns& mspc, const Double& time,
1476 : const Double& interval, const Int& antid) {
1477 : //Int start=lastIndex_p;
1478 176513 : Int start=lastIndexPerAnt_p[antid];
1479 176513 : Double tol=interval < 0.0 ? time*C::dbl_epsilon : interval/2.0;
1480 : // Search forwards
1481 176513 : Int nrows=mspc.time().nrow();
1482 956989 : for (Int i=start;i<nrows;i++) {
1483 : // Check for ANTENNA_ID. When antid<0 (default) ANTENNA_ID in POINTING table < 0 is ignored.
1484 : // MS v2 definition indicates ANTENNA_ID in POINING >= 0.
1485 956989 : if (antid >= 0 && mspc.antennaId()(i) != antid) {
1486 293484 : continue;
1487 : }
1488 :
1489 663505 : Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
1490 : // If the interval in the pointing table is negative, use the last
1491 : // entry. Note that this may be invalid (-1) but in that case
1492 : // the calling routine will generate an error
1493 663505 : Double mspc_interval = mspc.interval()(i);
1494 :
1495 663505 : if(mspc_interval<0.0) {
1496 : //return lastIndex_p;
1497 0 : return lastIndexPerAnt_p[antid];
1498 : }
1499 : // Pointing table interval is specified so we have to do a match
1500 : else {
1501 : // Is the midpoint of this pointing table entry within the specified
1502 : // tolerance of the main table entry?
1503 663505 : if(abs(midpoint-time) <= (mspc_interval/2.0+tol)) {
1504 : //lastIndex_p=i;
1505 176513 : lastIndexPerAnt_p[antid]=i;
1506 176513 : return i;
1507 : }
1508 : }
1509 : }
1510 : // Search backwards
1511 0 : for (Int i=start;i>=0;i--) {
1512 0 : if (antid >= 0 && mspc.antennaId()(i) != antid) {
1513 0 : continue;
1514 : }
1515 0 : Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
1516 0 : Double mspc_interval = mspc.interval()(i);
1517 0 : if(mspc_interval<0.0) {
1518 : //return lastIndex_p;
1519 0 : return lastIndexPerAnt_p[antid];
1520 : }
1521 : // Pointing table interval is specified so we have to do a match
1522 : else {
1523 : // Is the midpoint of this pointing table entry within the specified
1524 : // tolerance of the main table entry?
1525 0 : if (abs(midpoint-time) <= (mspc_interval/2.0+tol)) {
1526 : //lastIndex_p=i;
1527 0 : lastIndexPerAnt_p[antid]=i;
1528 0 : return i;
1529 : }
1530 : }
1531 : }
1532 : // No match!
1533 0 : return -1;
1534 : }
1535 :
1536 176495 : Bool SDGrid::getXYPos(const vi::VisBuffer2& vb, Int row) {
1537 :
1538 : Bool dointerp;
1539 176495 : const MSPointingColumns& act_mspc = vb.subtableColumns().pointing();
1540 176495 : Bool nullPointingTable = (act_mspc.nrow() < 1);
1541 176495 : Int pointIndex = -1;
1542 176495 : if (!nullPointingTable) {
1543 : ///if(vb.newMS()) vb.newMS does not work well using msid
1544 176495 : if (vb.msId() != msId_p) {
1545 139 : lastIndex_p = 0;
1546 139 : if (lastIndexPerAnt_p.nelements() < (size_t)vb.nAntennas()) {
1547 107 : lastIndexPerAnt_p.resize(vb.nAntennas());
1548 : }
1549 139 : lastIndexPerAnt_p = 0;
1550 139 : msId_p = vb.msId();
1551 139 : lastAntID_p = -1;
1552 : }
1553 176495 : pointIndex = getIndex(act_mspc, vb.time()(row), -1.0, vb.antenna1()(row));
1554 : //Try again to locate a pointing within the integration
1555 176495 : if (pointIndex < 0)
1556 0 : pointIndex = getIndex(act_mspc, vb.time()(row), vb.timeInterval()(row), vb.antenna1()(row));
1557 : }
1558 176495 : if (!nullPointingTable && ((pointIndex < 0) || (pointIndex >= Int(act_mspc.time().nrow())))) {
1559 0 : ostringstream o;
1560 0 : o << "Failed to find pointing information for time " <<
1561 0 : MVTime(vb.time()(row)/86400.0) << ": Omitting this point";
1562 0 : logIO_p << LogIO::DEBUGGING << String(o) << LogIO::POST;
1563 : // logIO_p << String(o) << LogIO::POST;
1564 0 : return false;
1565 : }
1566 :
1567 176495 : dointerp = false;
1568 176495 : if (!nullPointingTable && (vb.timeInterval()(row) < act_mspc.interval()(pointIndex))) {
1569 19602 : dointerp = true;
1570 19602 : if (!isSplineInterpolationReady) {
1571 16 : interpolator = new SDPosInterpolator(vb, pointingDirCol_p);
1572 16 : isSplineInterpolationReady = true;
1573 : } else {
1574 19586 : if (!interpolator->inTimeRange(vb.time()(row), vb.antenna1()(row))) {
1575 : // setup spline interpolator for the current dataset (CAS-11261, 2018/6/13 WK)
1576 188 : delete interpolator;
1577 188 : interpolator = 0;
1578 188 : interpolator = new SDPosInterpolator(vb, pointingDirCol_p);
1579 : }
1580 : }
1581 : }
1582 :
1583 176495 : if (!pointingToImage) {
1584 : // Set the frame
1585 124 : lastAntID_p = vb.antenna1()(row);
1586 248 : MPosition pos = vb.subtableColumns().antenna().positionMeas()(lastAntID_p);
1587 372 : MEpoch dummyEpoch(Quantity(0, "s"));
1588 124 : (!mFrame_p.epoch()) ? mFrame_p.set(dummyEpoch) : mFrame_p.resetEpoch(dummyEpoch);
1589 124 : (!mFrame_p.position()) ? mFrame_p.set(pos) : mFrame_p.resetPosition(pos);
1590 124 : if (!nullPointingTable) {
1591 124 : if (dointerp) {
1592 16 : worldPosMeas = directionMeas(act_mspc, pointIndex, vb.time()(row));
1593 : } else {
1594 108 : worldPosMeas = directionMeas(act_mspc, pointIndex);
1595 : }
1596 : } else {
1597 0 : worldPosMeas = vb.direction1()(row);
1598 : }
1599 :
1600 : // Make a machine to convert from the worldPosMeas to the output
1601 : // Direction Measure type for the relevant frame
1602 248 : MDirection::Ref outRef(directionCoord.directionType(), mFrame_p);
1603 124 : pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
1604 124 : if (!pointingToImage) {
1605 0 : logIO_p << "Cannot make direction conversion machine" << LogIO::EXCEPTION;
1606 : }
1607 :
1608 : // perform direction conversion to clear cache
1609 124 : MDirection _dir_tmp = (*pointingToImage)();
1610 : }
1611 :
1612 529485 : MEpoch epoch(Quantity(vb.time()(row), "s"));
1613 176495 : mFrame_p.resetEpoch(epoch);
1614 176495 : if (lastAntID_p != vb.antenna1()(row)) {
1615 14580 : if (lastAntID_p == -1) {
1616 : // antenna ID is unset
1617 15 : logIO_p << LogIO::DEBUGGING
1618 : << "update antenna position for conversion: new MS ID " << msId_p
1619 15 : << ", antenna ID " << vb.antenna1()(row) << LogIO::POST;
1620 : } else {
1621 14565 : logIO_p << LogIO::DEBUGGING
1622 : << "update antenna position for conversion: MS ID " << msId_p
1623 : << ", last antenna ID " << lastAntID_p
1624 14565 : << ", new antenna ID " << vb.antenna1()(row) << LogIO::POST;
1625 : }
1626 29160 : MPosition pos;
1627 14580 : lastAntID_p = vb.antenna1()(row);
1628 14580 : pos = vb.subtableColumns().antenna().positionMeas()(lastAntID_p);
1629 14580 : mFrame_p.resetPosition(pos);
1630 : }
1631 :
1632 176495 : if (!nullPointingTable) {
1633 176495 : if (dointerp) {
1634 39204 : MDirection newdir = directionMeas(act_mspc, pointIndex, vb.time()(row));
1635 : //Vector<Double> newdirv = newdir.getAngle("rad").getValue();
1636 19602 : worldPosMeas = (*pointingToImage)(newdir);
1637 : //cerr<<"dir0="<<newdirv(0)<<endl;
1638 :
1639 : //fprintf(pfile,"%.8f %.8f \n", newdirv(0), newdirv(1));
1640 : //printf("%lf %lf \n", newdirv(0), newdirv(1));
1641 : } else {
1642 156893 : worldPosMeas = (*pointingToImage)(directionMeas(act_mspc, pointIndex));
1643 : }
1644 : } else {
1645 0 : worldPosMeas = (*pointingToImage)(vb.direction1()(row));
1646 : }
1647 :
1648 176495 : Bool result = directionCoord.toPixel(xyPos, worldPosMeas);
1649 176495 : if (!result) {
1650 0 : logIO_p << "Failed to find a pixel for pointing direction of "
1651 0 : << MVTime(worldPosMeas.getValue().getLong("rad")).string(MVTime::TIME) << ", " << MVAngle(worldPosMeas.getValue().getLat("rad")).string(MVAngle::ANGLE) << LogIO::WARN << LogIO::POST;
1652 0 : return false;
1653 : }
1654 :
1655 176495 : if ((pointingDirCol_p == "SOURCE_OFFSET") || (pointingDirCol_p == "POINTING_OFFSET")) {
1656 : //there is no sense to track in offset coordinates...hopefully the
1657 : //user set the image coords right
1658 0 : fixMovingSource_p = false;
1659 : }
1660 176495 : if (fixMovingSource_p) {
1661 21292 : if (xyPosMovingOrig_p.nelements() < 2) {
1662 5 : directionCoord.toPixel(xyPosMovingOrig_p, firstMovingDir_p);
1663 : }
1664 : //via HADEC or AZEL for parallax of near sources
1665 42584 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
1666 42584 : MDirection tmphadec;
1667 21292 : if (upcase(movingDir_p.getRefString()).contains("APP")) {
1668 5828 : tmphadec = MDirection::Convert((vbutil_p->getEphemDir(vb, phaseCenterTime_p)), outref1)();
1669 15464 : } else if (upcase(movingDir_p.getRefString()).contains("COMET")) {
1670 2914 : MeasComet mcomet(Path(ephemTableName_p).absoluteName());
1671 2914 : mFrame_p.set(mcomet);
1672 2914 : tmphadec = MDirection::Convert(movingDir_p, outref1)();
1673 : } else {
1674 12550 : tmphadec = MDirection::Convert(movingDir_p, outref1)();
1675 : }
1676 42584 : MDirection actSourceDir = (*pointingToImage)(tmphadec);
1677 21292 : Vector<Double> actPix;
1678 21292 : directionCoord.toPixel(actPix, actSourceDir);
1679 :
1680 : //cout << row << " scan " << vb.scan()(row) << "xyPos " << xyPos << " xyposmovorig " << xyPosMovingOrig_p << " actPix " << actPix << endl;
1681 :
1682 21292 : xyPos = xyPos + xyPosMovingOrig_p - actPix;
1683 : }
1684 :
1685 176495 : return result;
1686 : // Convert to pixel coordinates
1687 : }
1688 :
1689 157037 : MDirection SDGrid::directionMeas(const MSPointingColumns& mspc, const Int& index){
1690 157037 : if (pointingDirCol_p == "TARGET") {
1691 0 : return mspc.targetMeas(index);
1692 157037 : } else if (pointingDirCol_p == "POINTING_OFFSET") {
1693 0 : if (!mspc.pointingOffsetMeasCol().isNull()) {
1694 0 : return mspc.pointingOffsetMeas(index);
1695 : }
1696 0 : cerr << "No PONTING_OFFSET column in POINTING table" << endl;
1697 157037 : } else if (pointingDirCol_p == "SOURCE_OFFSET") {
1698 0 : if (!mspc.sourceOffsetMeasCol().isNull()) {
1699 0 : return mspc.sourceOffsetMeas(index);
1700 : }
1701 0 : cerr << "No SOURCE_OFFSET column in POINTING table" << endl;
1702 157037 : } else if (pointingDirCol_p == "ENCODER") {
1703 0 : if (!mspc.encoderMeas().isNull()) {
1704 0 : return mspc.encoderMeas()(index);
1705 : }
1706 0 : cerr << "No ENCODER column in POINTING table" << endl;
1707 : }
1708 :
1709 : //default return this
1710 157037 : return mspc.directionMeas(index);
1711 : }
1712 :
1713 : // for the cases, interpolation of the pointing direction requires
1714 : // when data sampling rate higher than the pointing data recording
1715 : // (e.g. fast OTF)
1716 19618 : MDirection SDGrid::directionMeas(const MSPointingColumns& mspc, const Int& index,
1717 : const Double& time){
1718 : //spline interpolation
1719 19618 : if (isSplineInterpolationReady) {
1720 19618 : Int antid = mspc.antennaId()(index);
1721 19618 : if (interpolator->doSplineInterpolation(antid)) {
1722 19618 : return interpolator->interpolateDirectionMeasSpline(mspc, time, index, antid);
1723 : }
1724 : }
1725 :
1726 : //linear interpolation (as done before CAS-7911)
1727 : Int index1, index2;
1728 0 : if (time < mspc.time()(index)) {
1729 0 : if (index > 0) {
1730 0 : index1 = index-1;
1731 0 : index2 = index;
1732 0 : } else if (index == 0) {
1733 0 : index1 = index;
1734 0 : index2 = index+1;
1735 : }
1736 : } else {
1737 0 : if (index < Int(mspc.nrow()-1)) {
1738 0 : index1 = index;
1739 0 : index2 = index+1;
1740 0 : } else if (index == Int(mspc.nrow()-1) || (mspc.time()(index)-mspc.time()(index+1))>2*mspc.interval()(index)) {
1741 0 : index1 = index-1;
1742 0 : index2 = index;
1743 : }
1744 : }
1745 0 : return interpolateDirectionMeas(mspc, time, index, index1, index2);
1746 : }
1747 :
1748 0 : MDirection SDGrid::interpolateDirectionMeas(const MSPointingColumns& mspc,
1749 : const Double& time,
1750 : const Int& index,
1751 : const Int& indx1,
1752 : const Int& indx2){
1753 0 : Vector<Double> dir1,dir2;
1754 0 : Vector<Double> newdir(2),scanRate(2);
1755 : Double dLon, dLat;
1756 : Double ftime,ftime2,ftime1,dtime;
1757 0 : MDirection newDirMeas;
1758 0 : MDirection::Ref rf;
1759 : Bool isfirstRefPt;
1760 :
1761 0 : if (indx1 == index) {
1762 0 : isfirstRefPt = true;
1763 : } else {
1764 0 : isfirstRefPt = false;
1765 : }
1766 :
1767 0 : if (pointingDirCol_p == "TARGET") {
1768 0 : dir1 = mspc.targetMeas(indx1).getAngle("rad").getValue();
1769 0 : dir2 = mspc.targetMeas(indx2).getAngle("rad").getValue();
1770 0 : } else if (pointingDirCol_p == "POINTING_OFFSET") {
1771 0 : if (!mspc.pointingOffsetMeasCol().isNull()) {
1772 0 : dir1 = mspc.pointingOffsetMeas(indx1).getAngle("rad").getValue();
1773 0 : dir2 = mspc.pointingOffsetMeas(indx2).getAngle("rad").getValue();
1774 : } else {
1775 0 : cerr << "No PONTING_OFFSET column in POINTING table" << endl;
1776 : }
1777 0 : } else if (pointingDirCol_p == "SOURCE_OFFSET") {
1778 0 : if (!mspc.sourceOffsetMeasCol().isNull()) {
1779 0 : dir1 = mspc.sourceOffsetMeas(indx1).getAngle("rad").getValue();
1780 0 : dir2 = mspc.sourceOffsetMeas(indx2).getAngle("rad").getValue();
1781 : } else {
1782 0 : cerr << "No SOURCE_OFFSET column in POINTING table" << endl;
1783 : }
1784 0 : } else if (pointingDirCol_p == "ENCODER") {
1785 0 : if (!mspc.encoderMeas().isNull()) {
1786 0 : dir1 = mspc.encoderMeas()(indx1).getAngle("rad").getValue();
1787 0 : dir2 = mspc.encoderMeas()(indx2).getAngle("rad").getValue();
1788 : } else {
1789 0 : cerr << "No ENCODER column in POINTING table" << endl;
1790 : }
1791 : } else {
1792 0 : dir1 = mspc.directionMeas(indx1).getAngle("rad").getValue();
1793 0 : dir2 = mspc.directionMeas(indx2).getAngle("rad").getValue();
1794 : }
1795 :
1796 0 : dLon = dir2(0) - dir1(0);
1797 0 : dLat = dir2(1) - dir1(1);
1798 0 : ftime = floor(mspc.time()(indx1));
1799 0 : ftime2 = mspc.time()(indx2) - ftime;
1800 0 : ftime1 = mspc.time()(indx1) - ftime;
1801 0 : dtime = ftime2 - ftime1;
1802 0 : scanRate(0) = dLon/dtime;
1803 0 : scanRate(1) = dLat/dtime;
1804 : //scanRate(0) = dir2(0)/dtime-dir1(0)/dtime;
1805 : //scanRate(1) = dir2(1)/dtime-dir1(1)/dtime;
1806 : //Double delT = mspc.time()(index)-time;
1807 : //cerr<<"index="<<index<<" dLat="<<dLat<<" dtime="<<dtime<<" delT="<< delT<<endl;
1808 : //cerr<<"deldirlat="<<scanRate(1)*fabs(delT)<<endl;
1809 0 : if (isfirstRefPt) {
1810 0 : newdir(0) = dir1(0)+scanRate(0)*fabs(mspc.time()(index)-time);
1811 0 : newdir(1) = dir1(1)+scanRate(1)*fabs(mspc.time()(index)-time);
1812 0 : rf = mspc.directionMeas(indx1).getRef();
1813 : } else {
1814 0 : newdir(0) = dir2(0)-scanRate(0)*fabs(mspc.time()(index)-time);
1815 0 : newdir(1) = dir2(1)-scanRate(1)*fabs(mspc.time()(index)-time);
1816 0 : rf = mspc.directionMeas(indx2).getRef();
1817 : }
1818 : //default return this
1819 0 : Quantity rDirLon(newdir(0), "rad");
1820 0 : Quantity rDirLat(newdir(1), "rad");
1821 0 : newDirMeas = MDirection(rDirLon, rDirLat, rf);
1822 : //cerr<<"newDirMeas rf="<<newDirMeas.getRefString()<<endl;
1823 : //return mspc.directionMeas(index);
1824 0 : return newDirMeas;
1825 : }
1826 :
1827 169210 : void SDGrid::pickWeights(const vi::VisBuffer2& vb, Matrix<Float>& weight){
1828 : //break reference
1829 169210 : weight.resize();
1830 :
1831 169210 : if (useImagingWeight_p) {
1832 0 : weight.reference(vb.imagingWeight());
1833 : } else {
1834 338420 : const Cube<Float> weightspec(vb.weightSpectrum());
1835 169210 : weight.resize(vb.nChannels(), vb.nRows());
1836 :
1837 169210 : if (weightspec.nelements() == 0) {
1838 0 : for (rownr_t k = 0; k < vb.nRows(); ++k) {
1839 : //cerr << "nrow " << vb.nRow() << " " << weight.shape() << " " << weight.column(k).shape() << endl;
1840 0 : weight.column(k).set(mean(vb.weight().column(k)));
1841 : }
1842 : } else {
1843 169210 : Int npol = weightspec.shape()(0);
1844 169210 : if (npol == 1) {
1845 100 : for (rownr_t k = 0; k < vb.nRows(); ++k) {
1846 106 : for (int chan = 0; chan < vb.nChannels(); ++chan) {
1847 56 : weight(chan, k)=weightspec(0, chan, k);
1848 : }
1849 : }
1850 : } else {
1851 345605 : for (rownr_t k = 0; k < vb.nRows(); ++k) {
1852 94810341 : for (int chan = 0; chan < vb.nChannels(); ++chan) {
1853 94633896 : weight(chan, k) = (weightspec(0, chan, k) + weightspec((npol-1), chan, k))/2.0f;
1854 : }
1855 : }
1856 : }
1857 : }
1858 : }
1859 169210 : }
1860 :
1861 129 : void SDGrid::clipMinMax() {
1862 129 : if (clipminmax_) {
1863 : Bool gmin_b, gmax_b, wmin_b, wmax_b, np_b;
1864 8 : const auto *gmin_p = gmin_.getStorage(gmin_b);
1865 8 : const auto *gmax_p = gmax_.getStorage(gmax_b);
1866 8 : const auto *wmin_p = wmin_.getStorage(wmin_b);
1867 8 : const auto *wmax_p = wmax_.getStorage(wmax_b);
1868 8 : const auto *np_p = npoints_.getStorage(np_b);
1869 :
1870 : Bool data_b, weight_b, sumw_b;
1871 8 : auto data_p = griddedData.getStorage(data_b);
1872 8 : auto weight_p = wGriddedData.getStorage(weight_b);
1873 8 : auto sumw_p = sumWeight.getStorage(sumw_b);
1874 :
1875 16 : auto arrayShape = griddedData.shape();
1876 8 : size_t num_xy = arrayShape.getFirst(2).product();
1877 8 : size_t num_polchan = arrayShape.getLast(2).product();
1878 80 : for (size_t i = 0; i < num_xy; ++i) {
1879 153 : for (size_t j = 0; j < num_polchan; ++j) {
1880 81 : auto k = i * num_polchan + j;
1881 81 : if (np_p[k] == 1) {
1882 2 : auto wt = wmin_p[k];
1883 2 : data_p[k] = wt * gmin_p[k];
1884 2 : weight_p[k] = wt;
1885 2 : sumw_p[j] += wt;
1886 79 : } else if (np_p[k] == 2) {
1887 1 : auto wt = wmin_p[k];
1888 1 : data_p[k] = wt * gmin_p[k];
1889 1 : weight_p[k] = wt;
1890 1 : sumw_p[j] += wt;
1891 1 : wt = wmax_p[k];
1892 1 : data_p[k] += wt * gmax_p[k];
1893 1 : weight_p[k] += wt;
1894 1 : sumw_p[j] += wt;
1895 : }
1896 : }
1897 : }
1898 :
1899 8 : wGriddedData.putStorage(weight_p, weight_b);
1900 8 : griddedData.putStorage(data_p, data_b);
1901 8 : sumWeight.putStorage(sumw_p, sumw_b);
1902 :
1903 8 : npoints_.freeStorage(np_p, np_b);
1904 8 : wmax_.freeStorage(wmax_p, wmax_b);
1905 8 : wmin_.freeStorage(wmin_p, wmin_b);
1906 8 : gmax_.freeStorage(gmax_p, gmax_b);
1907 8 : gmin_.freeStorage(gmin_p, gmin_b);
1908 : }
1909 129 : }
1910 :
1911 : } //End of namespace refim
1912 : } //#End casa namespace
|
