Line data Source code
1 : //# SimplePBConvFunc.cc: implementation of SimplePBConvFunc
2 : //# Copyright (C) 2007-2020
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 adressed 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 : //#
27 : //# $Id$
28 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <casacore/casa/Arrays/Array.h>
30 : #include <casacore/casa/Arrays/MaskedArray.h>
31 : #include <casacore/casa/Arrays/Vector.h>
32 : #include <casacore/casa/Arrays/Slice.h>
33 : #include <casacore/casa/Arrays/Slicer.h>
34 : #include <casacore/casa/Arrays/Matrix.h>
35 : #include <casacore/casa/Arrays/Cube.h>
36 : #include <casacore/casa/OS/Timer.h>
37 : #include <casacore/casa/Utilities/Assert.h>
38 : #include <casacore/casa/Quanta/MVTime.h>
39 : #include <casacore/casa/Quanta/MVAngle.h>
40 : #include <casacore/measures/Measures/MeasTable.h>
41 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
42 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
43 :
44 : #include <casacore/images/Images/ImageInterface.h>
45 : #include <casacore/images/Images/PagedImage.h>
46 : #include <casacore/images/Images/TempImage.h>
47 : #include <casacore/images/Images/SubImage.h>
48 : #include <casacore/casa/Logging/LogIO.h>
49 : #include <casacore/casa/Logging/LogSink.h>
50 : #include <casacore/casa/Logging/LogMessage.h>
51 :
52 : #include <casacore/ms/MeasurementSets/MSColumns.h>
53 : #include <casacore/lattices/Lattices/ArrayLattice.h>
54 : #include <casacore/lattices/Lattices/SubLattice.h>
55 : #include <casacore/lattices/LRegions/LCBox.h>
56 : #include <casacore/lattices/LEL/LatticeExpr.h>
57 : #include <casacore/lattices/Lattices/LatticeCache.h>
58 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
59 :
60 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
61 :
62 : #include <msvis/MSVis/VisBuffer.h>
63 : #include <msvis/MSVis/VisibilityIterator.h>
64 :
65 : #include <synthesis/TransformMachines2/SimplePBConvFunc.h>
66 : #include <synthesis/TransformMachines2/SkyJones.h>
67 : #include <synthesis/TransformMachines2/AWConvFuncHolder.h>
68 : #include <casacore/casa/Utilities/CompositeNumber.h>
69 : #include <iomanip>
70 : #include <math.h>
71 :
72 : using namespace casacore;
73 : namespace casa { //# NAMESPACE CASA - BEGIN
74 : namespace refim {//# namespace for imaging refactor
75 : using namespace casacore;
76 : using namespace casa;
77 : using namespace casacore;
78 : using namespace casa::refim;
79 :
80 0 : SimplePBConvFunc::SimplePBConvFunc() : nchan_p(-1),
81 0 : npol_p(-1), pointToPix_p(), directionIndex_p(-1), thePix_p(0),
82 0 : filledFluxScale_p(false), doneMainConv_p(0),
83 :
84 0 : calcFluxScale_p(true), usePointingTable_p(False), actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), pointingPix_p(), awConvs_p(nullptr) {
85 : //
86 :
87 0 : pbClass_p = PBMathInterface::COMMONPB;
88 0 : ft_p = FFT2D(true);
89 0 : usePointingTable_p = False;
90 0 : }
91 :
92 0 : SimplePBConvFunc::SimplePBConvFunc(const PBMathInterface::PBClass typeToUse) :
93 0 : nchan_p(-1), npol_p(-1), pointToPix_p(),
94 0 : directionIndex_p(-1), thePix_p(0), filledFluxScale_p(false), doneMainConv_p(0),
95 0 : calcFluxScale_p(true), usePointingTable_p(False), actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), pointingPix_p(), awConvs_p(nullptr) {
96 : //
97 0 : pbClass_p = typeToUse;
98 0 : ft_p = FFT2D(true);
99 0 : usePointingTable_p = False;
100 0 : }
101 :
102 0 : SimplePBConvFunc::SimplePBConvFunc(const RecordInterface& rec, const Bool calcfluxneeded)
103 0 : : nchan_p(-1), npol_p(-1), pointToPix_p(), directionIndex_p(-1), thePix_p(0), filledFluxScale_p(false),
104 0 : doneMainConv_p(0),
105 0 : calcFluxScale_p(calcfluxneeded), usePointingTable_p(False), actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), pointingPix_p(), awConvs_p(nullptr)
106 : {
107 0 : String err;
108 0 : fromRecord(err, rec, calcfluxneeded);
109 0 : ft_p = FFT2D(true);
110 0 : usePointingTable_p = False;
111 0 : }
112 :
113 0 : SimplePBConvFunc::~SimplePBConvFunc() {
114 : //
115 :
116 0 : }
117 :
118 0 : void SimplePBConvFunc::storeImageParams(const ImageInterface<Complex>& iimage,
119 : const vi::VisBuffer2& vb) {
120 : //image signature changed...rather simplistic for now
121 0 : if ((iimage.shape().product() != nx_p * ny_p * nchan_p * npol_p) || nchan_p < 1) {
122 0 : csys_p = iimage.coordinates();
123 0 : Int coordIndex = csys_p.findCoordinate(Coordinate::DIRECTION);
124 0 : AlwaysAssert(coordIndex >= 0, AipsError);
125 0 : directionIndex_p = coordIndex;
126 0 : dc_p = csys_p.directionCoordinate(directionIndex_p);
127 0 : ObsInfo imInfo = csys_p.obsInfo();
128 0 : String tel = imInfo.telescope();
129 0 : MPosition pos;
130 0 : MSColumns mscol(vb.ms());
131 0 : if (vb.subtableColumns().observation().nrow() > 0) {
132 0 : tel = vb.subtableColumns().observation().telescopeName()(mscol.observationId()(0));
133 : }
134 0 : if (tel.length() == 0 || !tel.contains("VLA") ||
135 0 : !MeasTable::Observatory(pos, tel)) {
136 : // unknown observatory, use first antenna
137 0 : Int ant1 = vb.antenna1()(0);
138 0 : pos = vb.subtableColumns().antenna().positionMeas()(ant1);
139 : }
140 0 : imInfo.setTelescope(tel);
141 0 : csys_p.setObsInfo(imInfo);
142 : //Store this to build epochs via the time access of visbuffer later
143 0 : timeMType_p = MEpoch::castType(mscol.timeMeas()(0).getRef().getType());
144 0 : timeUnit_p = Unit(mscol.timeMeas().measDesc().getUnits()(0).getName());
145 : // timeUnit_p=Unit("s");
146 : //cout << "UNIT " << timeUnit_p.getValue() << " name " << timeUnit_p.getName() << endl;
147 0 : pointFrame_p = MeasFrame(imInfo.obsDate(), pos);
148 0 : MDirection::Ref elRef(dc_p.directionType(), pointFrame_p);
149 : //For now we set the conversion from this direction
150 0 : pointToPix_p = MDirection::Convert(MDirection(), elRef);
151 0 : nx_p = iimage.shape()(coordIndex);
152 0 : ny_p = iimage.shape()(coordIndex + 1);
153 0 : pointingPix_p.resize(nx_p, ny_p);
154 0 : pointingPix_p.set(false);
155 0 : coordIndex = csys_p.findCoordinate(Coordinate::SPECTRAL);
156 0 : Int pixAxis = csys_p.pixelAxes(coordIndex)[0];
157 0 : nchan_p = iimage.shape()(pixAxis);
158 0 : coordIndex = csys_p.findCoordinate(Coordinate::STOKES);
159 0 : pixAxis = csys_p.pixelAxes(coordIndex)[0];
160 0 : npol_p = iimage.shape()(pixAxis);
161 0 : if (calcFluxScale_p) {
162 0 : if (fluxScale_p.shape().nelements() == 0) {
163 0 : fluxScale_p = TempImage<Float>(IPosition(4, nx_p, ny_p, npol_p, nchan_p), csys_p);
164 0 : fluxScale_p.set(0.0);
165 : }
166 0 : filledFluxScale_p = false;
167 : }
168 :
169 0 : }
170 :
171 0 : }
172 :
173 0 : void SimplePBConvFunc::toPix(const vi::VisBuffer2& vb, const MVDirection& extraShift, const Bool useExtraShift) {
174 0 : thePix_p.resize(2);
175 :
176 0 : const MDirection& p1 = pointingDirAnt1(vb);
177 0 : if (dc_p.directionType() != MDirection::castType(p1.getRef().getType())) {
178 : //pointToPix_p.setModel(theDir);
179 :
180 0 : String tel = csys_p.obsInfo().telescope();
181 0 : if (!tel.contains("VLA")) {
182 : //use first antenna as direction1_p is used to calculate pointing
183 : // as only VLA uses observatory pos for calculations
184 0 : Int ant1 = vb.antenna1()(0);
185 0 : MPosition pos = vb.subtableColumns().antenna().positionMeas()(ant1);
186 0 : pointFrame_p.resetPosition(pos);
187 0 : }
188 0 : MEpoch timenow(Quantity(vb.time()(0), timeUnit_p), timeMType_p);
189 : //cerr << "Ref " << vb.direction1()(0).getRefString() << " ep " << timenow.getRefString() << " time " << MVTime(timenow.getValue().getTime()).string(MVTime::YMD) << endl;
190 0 : pointFrame_p.resetEpoch(timenow);
191 : ///////////////////////////
192 : //MDirection some=pointToPix_p(vb.direction1()(0));
193 : //MVAngle mvRa=some.getAngle().getValue()(0);
194 : //MVAngle mvDec=some.getAngle().getValue()(1);
195 :
196 : //cout << mvRa(0.0).string(MVAngle::TIME,8) << " ";
197 : // cout << mvDec.string(MVAngle::DIG2,8) << " ";
198 : //cout << MDirection::showType(some.getRefPtr()->getType()) << endl;
199 :
200 : //////////////////////////
201 : //pointToPix holds pointFrame_p by reference...
202 : //thus good to go for conversion
203 0 : direction1_p = pointToPix_p(p1);
204 : //direction2_p=pointToPix_p(vb.direction2()(0));
205 0 : direction2_p = direction1_p;
206 :
207 :
208 0 : }
209 : else {
210 0 : direction1_p = p1;
211 :
212 : //direction2_p=vb.direction2()(0);
213 : //For now
214 0 : direction2_p = direction1_p;
215 :
216 : }
217 : //Should return both antennas direction in the future
218 :
219 0 : if (useExtraShift) {
220 0 : direction1_p.shift(extraShift, False);
221 0 : direction2_p.shift(extraShift, False);
222 : }
223 0 : dc_p.toPixel(thePix_p, direction1_p);
224 :
225 0 : }
226 :
227 0 : void SimplePBConvFunc::setWeightImage(CountedPtr<TempImage<Float> >& wgtimage) {
228 0 : convWeightImage_p = wgtimage;
229 0 : filledFluxScale_p = false;
230 0 : calcFluxScale_p = true;
231 :
232 0 : }
233 :
234 0 : void SimplePBConvFunc::reset() {
235 0 : doneMainConv_p.resize();
236 0 : convFunctions_p.resize(0, true);
237 0 : convWeights_p.resize(0, true);
238 0 : convSizes_p.resize(0, true);
239 0 : convSupportBlock_p.resize(0, true);
240 0 : convFunctionMap_p.clear();
241 0 : vbConvIndex_p.clear();
242 0 : awConvs_p=nullptr;
243 0 : ft_p = FFT2D(true);
244 0 : }
245 :
246 :
247 :
248 0 : Int SimplePBConvFunc::convIndex(const vi::VisBuffer2& vb, const uInt nchan) {
249 0 : String elkey = String::toString(vb.msId()) + String("_") + String::toString(vb.spectralWindows()[0]) + String("_") + String::toString(nchan);
250 0 : if (vbConvIndex_p.count(elkey) > 0) {
251 0 : return vbConvIndex_p[elkey];
252 : }
253 0 : Int val = vbConvIndex_p.size();
254 0 : vbConvIndex_p[elkey] = val;
255 0 : return val;
256 0 : }
257 :
258 0 : const MDirection& SimplePBConvFunc::pointingDirAnt1(const vi::VisBuffer2& vb) {
259 :
260 :
261 0 : std::ostringstream oss;
262 :
263 0 : oss << vb.msId() << "_" << vb.antenna1()(0) << "_";
264 0 : oss.precision(13);
265 0 : oss << vb.time()(0);
266 0 : String elkey = oss.str();
267 : // String elkey=String::toString(vb.msId())+String("_")+String::toString(vb.antenna1()(0))+String("_")
268 : // +String::toString(vb.time()(0));
269 :
270 : //cerr << "key " << elkey << " count " << ant1PointVal_p.count(elkey) << " size " << ant1PointVal_p.size() << " " << ant1PointingCache_p.nelements() << endl;
271 0 : if (ant1PointVal_p.count(elkey) > 0) {
272 0 : return ant1PointingCache_p[ant1PointVal_p[elkey]];
273 :
274 : }
275 0 : Bool hasValidPointing = False;
276 0 : if (Table::isReadable(vb.ms().pointingTableName())) {
277 0 : hasValidPointing = usePointingTable_p && (vb.ms().pointing().nrow() > 0);
278 : }
279 :
280 0 : Int val = ant1PointingCache_p.nelements();
281 0 : ant1PointingCache_p.resize(val + 1, true);
282 0 : if (hasValidPointing) {
283 : //ant1PointingCache_p[val]=vb.direction1()[0];
284 0 : ant1PointingCache_p[val] = vbutil_p->getPointingDir(vb, vb.antenna1()(0), 0, dc_p.directionType());
285 : }
286 : else
287 0 : ant1PointingCache_p[val] = vbutil_p->getPhaseCenter(vb);
288 : //ant1PointingCache_p[val]=vbUtil_p.getPointingDir(vb, vb.antenna1()(0), 0);
289 0 : ant1PointVal_p[elkey] = val;
290 0 : return ant1PointingCache_p[val];
291 :
292 0 : }
293 :
294 0 : void SimplePBConvFunc::rephaseConvFunc(const ImageInterface<Complex>& iimage,
295 : const vi::VisBuffer2& vb, const Int& convSampling, Array<Complex>& convFunc,
296 : Array<Complex>& weightConvFunc, const std::vector<casacore::Int>& pmap,
297 : const std::vector<casacore::Int>& cmap,
298 : const std::vector<casacore::Int>& rmap, const MVDirection& extraShift, const Bool useExtraShift) {
299 : /*storeImageParams(iimage,vb);
300 : toPix(vb, extraShift, useExtraShift);
301 : Vector<Double> pixFieldDir(2);
302 : pixFieldDir=thePix_p;
303 : pixFieldDir(0)=pixFieldDir(0)- Double(nx / 2);
304 : pixFieldDir(1)=pixFieldDir(1)- Double(ny / 2);
305 : pixFieldDir(0)=-pixFieldDir(0)*2.0*M_PI/Double(nx)/Double(convSamp)ling;
306 : pixFieldDir(1)=-pixFieldDir(1)*2.0*M_PI/Double(ny)/Double(convSampling);
307 : Int nconvrow=convFunc.shape()(4);
308 : Int nconvchan=convFunc.shape(3);
309 : Int nconvpol=convFunc.shape()(2);
310 : Int convsize=convFunc.shape()(0);
311 : Bool delc;
312 : Bool delw;
313 : Double dirX=pixFieldDir(0);
314 : Double dirY=pixFieldDir(1);
315 : Complex *convstor=convFunc.getStorage(delc);
316 : Complex *weightstor=weightConvFunc.getStorage(delw);
317 : #pragma omp parallel default(none) firstprivate(convstor, weightstor, dirX, dirY, convsize, nconvrow, nconvchan, nconvpol)
318 : {
319 :
320 : #pragma omp for
321 : for(Int iy=0; iy<convsize; ++iy) {
322 : applyGradientToYLine(iy, convstor, weightstor, dirX, dirY, convsize, nconvrow, nconvchan, nconvpol);
323 :
324 : }
325 : }///End of pragma
326 : convFunc.putStorage(convstor, delc);
327 : weightConvFunc.putStorage(weightstor, delw);
328 : */
329 0 : throw(AipsError("Programmers' error: no implemented"));
330 : }
331 0 : void SimplePBConvFunc::findConvFunction(const ImageInterface<Complex>& iimage,
332 : const vi::VisBuffer2& vb,
333 : const Int& convSampling,
334 : const Vector<Double>& visFreq,
335 : Array<Complex>& convFunc,
336 : Array<Complex>& weightConvFunc,
337 : Vector<Int>& convsize,
338 : Vector<Int>& convSupport,
339 : Vector<Int>& convFuncPolMap,
340 : Vector<Int>& convFuncChanMap,
341 : Vector<Int>& convFuncRowMap,
342 : const Bool /*getConjFreqConvFunc*/,
343 : const MVDirection& extraShift, const Bool useExtraShift
344 : ) {
345 :
346 :
347 :
348 0 : Int convSamp = 2 * convSampling;
349 : /////////////////////////
350 0 : LogIO os1;
351 0 : os1 << "msID " << vb.msId() << " ANT1 id" << vb.antenna1()(0) << " direction " << vb.direction1()(0).toString() << " ANT2 id" << vb.antenna2()(0) << " direction " << vb.direction2()(0).toString() << LogIO::DEBUG1 << LogIO::POST;
352 : //////////////////////
353 0 : storeImageParams(iimage, vb);
354 0 : convFuncChanMap.resize(vb.nChannels());
355 0 : Vector<Double> beamFreqs;
356 0 : findUsefulChannels(convFuncChanMap, beamFreqs, vb, visFreq);
357 : //cerr << "CHANMAP " << convFuncChanMap << " beamFreqs " << beamFreqs << endl;
358 0 : Int nBeamChans = beamFreqs.nelements();
359 : //indgen(convFuncChanMap);
360 0 : convFuncPolMap.resize(vb.nCorrelations());
361 0 : convFuncPolMap.set(0);
362 : //Only one plane in this version
363 0 : convFuncRowMap.resize();
364 0 : convFuncRowMap = Vector<Int>(vb.nRows(), 0);
365 : //break reference
366 0 : convFunc.resize();
367 0 : weightConvFunc.resize();
368 0 : LogIO os;
369 0 : os << LogOrigin("SimplePBConv", "findConvFunction") << LogIO::NORMAL;
370 :
371 :
372 : // Get the coordinate system
373 0 : CoordinateSystem coords(iimage.coordinates());
374 :
375 :
376 0 : actualConvIndex_p = convIndex(vb, visFreq.nelements());
377 : //cerr << "In findConv " << actualConvIndex_p << endl;
378 : // Make a two dimensional image to calculate the
379 : // primary beam. We want this on a fine grid in the
380 : // UV plane
381 0 : Int directionIndex = directionIndex_p;
382 0 : AlwaysAssert(directionIndex >= 0, AipsError);
383 :
384 : // Set up the convolution function.
385 0 : Int nx = nx_p;
386 0 : Int ny = ny_p;
387 : // convSize_p=max(nx,ny)*convSampling;
388 : //cerr << "size " << nx << " " << ny << endl;
389 : //3 times the support size
390 0 : if (doneMainConv_p.shape()[0] < (actualConvIndex_p + 1)) {
391 : // cerr << "resizing DONEMAIN " << doneMainConv_p.shape()[0] << endl;
392 0 : doneMainConv_p.resize(actualConvIndex_p + 1, true);
393 0 : doneMainConv_p[actualConvIndex_p] = false;
394 : }
395 :
396 0 : if (!(doneMainConv_p[actualConvIndex_p])) {
397 :
398 : //convSize_p=4*(sj_p->support(vb, coords));
399 0 : convSize_p = Int(max(nx_p, ny_p) / 2) * 2 * convSamp;
400 : // Make this a nice composite number, to speed up FFTs
401 : //cerr << "convSize_p 0 " << convSize_p << " convSamp " << convSamp<< endl;
402 0 : CompositeNumber cn(uInt(convSize_p * 2.0));
403 :
404 0 : convSize_p = cn.nextLargerEven(Int(convSize_p));
405 : //cerr << "convSize : " << convSize_p << endl;
406 :
407 0 : }
408 :
409 :
410 0 : toPix(vb, extraShift, useExtraShift);
411 : //Timer tim;
412 : //tim.mark();
413 0 : addPBToFlux(vb);
414 : //tim.show("After addPBToFlux");
415 0 : DirectionCoordinate dc = dc_p;
416 :
417 : //where in the image in pixels is this pointing
418 0 : Vector<Double> pixFieldDir(2);
419 0 : pixFieldDir = thePix_p;
420 :
421 : //cerr << "pix of pointing " << pixFieldDir << endl;
422 0 : MDirection fieldDir = direction1_p;
423 : //shift from center
424 0 : pixFieldDir(0) = pixFieldDir(0) - Double(nx / 2);
425 0 : pixFieldDir(1) = pixFieldDir(1) - Double(ny / 2);
426 :
427 : //phase gradient per pixel to apply
428 0 : pixFieldDir(0) = -pixFieldDir(0) * 2.0 * M_PI / Double(nx) / Double(convSampling);
429 0 : pixFieldDir(1) = -pixFieldDir(1) * 2.0 * M_PI / Double(ny) / Double(convSampling);
430 :
431 : //cerr << "DonemainConv " << doneMainConv_p[actualConvIndex_p] << endl;
432 0 : if (!doneMainConv_p[actualConvIndex_p]) {
433 0 : Vector<Double> sampling;
434 0 : sampling = dc.increment();
435 0 : sampling *= Double(convSamp);
436 0 : sampling(0) *= Double(nx) / Double(convSize_p);
437 0 : sampling(1) *= Double(ny) / Double(convSize_p);
438 0 : dc.setIncrement(sampling);
439 :
440 :
441 0 : Vector<Double> unitVec(2);
442 0 : unitVec = convSize_p / 2;
443 0 : dc.setReferencePixel(unitVec);
444 :
445 :
446 : //make sure we are using the same units
447 0 : fieldDir.set(dc.worldAxisUnits()(0));
448 :
449 0 : dc.setReferenceValue(fieldDir.getAngle().getValue());
450 :
451 0 : coords.replaceCoordinate(dc, directionIndex);
452 0 : Int spind = coords.findCoordinate(Coordinate::SPECTRAL);
453 0 : SpectralCoordinate spCoord = coords.spectralCoordinate(spind);
454 0 : spCoord.setReferencePixel(Vector<Double>(1, 0.0));
455 0 : spCoord.setReferenceValue(Vector<Double>(1, beamFreqs(0)));
456 0 : if (beamFreqs.nelements() > 1)
457 0 : spCoord.setIncrement(Vector<Double>(1, beamFreqs(1) - beamFreqs(0)));
458 0 : coords.replaceCoordinate(spCoord, spind);
459 :
460 :
461 0 : CoordinateSystem coordLastPlane = coords;
462 0 : spCoord.setReferenceValue(Vector<Double>(1, beamFreqs(nBeamChans - 1)));
463 0 : coordLastPlane.replaceCoordinate(spCoord, spind);
464 : //cerr << "BEAM freqs " << beamFreqs << endl;
465 :
466 : // coords.list(logIO(), MDoppler::RADIO, IPosition(), IPosition());
467 :
468 0 : Int tempConvSize = ((convSize_p / 4 / (convSamp / convSampling)) / 2) * 2;
469 0 : IPosition pbShape(4, tempConvSize, tempConvSize, 1, nBeamChans);
470 :
471 0 : Long memtot = HostInfo::memoryFree();
472 0 : Double memtobeused = Double(memtot) * 1024.0;
473 : //check for 32 bit OS and limit it to 2Gbyte
474 : if (sizeof(void*) == 4) {
475 : if (memtot > 2000000)
476 : memtot = 2000000;
477 : }
478 0 : if (memtot <= 2000000)
479 0 : memtobeused = 0.0;
480 : //cerr << "mem to be used " << memtobeused << endl;
481 : //tim.mark();
482 0 : IPosition start(4, 0, 0, 0, 0);
483 : //IPosition pbSlice(4, convSize_p, convSize_p, 1, 1);
484 : //cerr << "pbshape " << pbShape << endl;
485 0 : TempImage<Complex> twoDPB(TiledShape(pbShape, IPosition(4, pbShape(0), pbShape(1), 1, 1)), coords, memtobeused / 10.0);
486 :
487 : //tim.show("after making one image");
488 0 : convFunc_p.resize(tempConvSize, tempConvSize);
489 0 : convFunc_p = 0.0;
490 :
491 :
492 :
493 : // Accumulate terms
494 : //Matrix<Complex> screen(convSize_p, convSize_p);
495 : //screen=1.0;
496 : // Either the SkyJones
497 : //tim.mark();
498 : //twoDPB.set(Complex(1.0,0.0));
499 : //for (Int k=0; k < nBeamChans; ++k){
500 : //blcin[3]=k;
501 : //trcin[3]=k;
502 : //Slicer slin(blcin, trcin, Slicer::endIsLast);
503 : //SubImage<Complex> subim(twoDPB, slin, true);
504 0 : TempImage<Complex> subim(IPosition(4, convSize_p, convSize_p, 1, 1), coordLastPlane, memtobeused / 2.2);
505 0 : subim.set(Complex(1.0, 0.0));
506 : //twoDPB.putSlice(screen, start);
507 0 : sj_p->apply(subim, subim, vb, 0);
508 : //LatticeFFT::cfft2d(subim);
509 0 : ft_p.c2cFFT(subim);
510 : // }
511 : //tim.show("after an apply" );
512 : //tim.mark();
513 0 : TempImage<Float> screen2(TiledShape(IPosition(4, convSize_p, convSize_p, 1, 1)), coordLastPlane, memtobeused / 2.2);
514 0 : screen2.set(1.0);
515 0 : TempImage<Complex> subout(TiledShape(IPosition(4, convSize_p, convSize_p, 1, 1)), coordLastPlane, memtobeused / 2.2);
516 : //////Making a reference on half of the lattice as on the Mac rcfft is failing for some
517 : //////reason
518 0 : SubImage<Complex> halfsubout(subout, Slicer(IPosition(4, 0, 0, 0, 0), IPosition(4, convSize_p / 2, convSize_p - 1, 0, 0), Slicer::endIsLast), true);
519 0 : sj_p->applySquare(screen2, screen2, vb, 0);
520 0 : ft_p.r2cFFT(halfsubout, screen2);
521 : //LatticeFFT::rcfft(halfsubout, screen2, true, false);
522 : //Real FFT fills only first half of the array
523 : //making it look like a Complex to Complex FFT
524 0 : IPosition iblc(4, 0, 3 * subout.shape()(1) / 8, 0, 0);
525 0 : IPosition itrc(4, 0, 5 * subout.shape()(1) / 8, 0, 0);
526 0 : for (Int x = subout.shape()(0) / 2; x < (5 * subout.shape()(0) / 8); ++x) {
527 :
528 0 : iblc[0] = x - subout.shape()(0) / 2;
529 0 : itrc[0] = x - subout.shape()(0) / 2;
530 0 : Slicer isl(iblc, itrc, Slicer::endIsLast);
531 0 : iblc[0] = x;
532 0 : subout.putSlice(subout.getSlice(isl), iblc);
533 0 : }
534 0 : for (Int x = subout.shape()(0) / 2 + 1; x < (5 * subout.shape()(0) / 8); ++x) {
535 :
536 0 : iblc[0] = x;
537 0 : itrc[0] = x;
538 0 : Slicer isl(iblc, itrc, Slicer::endIsLast);
539 0 : iblc[0] = subout.shape()(0) - x;
540 0 : subout.putSlice(subout.getSlice(isl), iblc);
541 0 : if (x == (subout.shape()(0) - 1)) {
542 0 : iblc[0] = 0;
543 0 : subout.putSlice(subout.getSlice(isl), iblc);
544 : }
545 0 : }
546 : //End of FFT's
547 : //tim.show("After apply2 ");
548 0 : TempImage<Complex> twoDPB2(TiledShape(pbShape, IPosition(4, pbShape(0), pbShape(1), 1, 1)), coords, memtobeused / 10.0);
549 :
550 0 : IPosition blcout(4, 0, 0, 0, nBeamChans - 1);
551 0 : IPosition trcout(4, pbShape(0) - 1, pbShape(1) - 1, 0, nBeamChans - 1);
552 0 : Slicer outsl(blcout, trcout, Slicer::endIsLast);
553 :
554 0 : IPosition blcin(4, convSize_p / 2 - pbShape(0) / 2, convSize_p / 2 - pbShape(1) / 2, 0, 0);
555 0 : IPosition trcin(4, convSize_p / 2 + pbShape(0) / 2 - 1, convSize_p / 2 + pbShape(1) / 2 - 1, 0, 0);
556 0 : Slicer insl(blcin, trcin, Slicer::endIsLast);
557 : {
558 0 : SubImage<Complex> subtwoDPB(twoDPB, outsl, true);
559 0 : SubImage<Complex> intwoDPB(subim, insl, false);
560 0 : subtwoDPB.copyData(intwoDPB);
561 0 : }
562 : {
563 0 : SubImage<Complex> subtwoDPB2(twoDPB2, outsl, true);
564 0 : SubImage<Complex> intwoDPB2(subout, insl, false);
565 0 : subtwoDPB2.copyData(intwoDPB2);
566 0 : }
567 :
568 0 : if (nBeamChans > 0) {
569 0 : blcin = IPosition(4, 0, 0, 0, nBeamChans - 1);
570 0 : trcin = IPosition(4, pbShape(0) - 1, pbShape(1) - 1, 0, nBeamChans - 1);
571 0 : Slicer slin(blcin, trcin, Slicer::endIsLast);
572 0 : SubImage<Complex> origPB(twoDPB, slin, false);
573 0 : IPosition elshape = origPB.shape();
574 0 : Matrix<Complex> i1 = origPB.get(true);
575 0 : SubImage<Complex> origPB2(twoDPB2, slin, false);
576 0 : Matrix<Complex> i2 = origPB2.get(true);
577 0 : Int cenX = i1.shape()(0) / 2;
578 0 : Int cenY = i1.shape()(1) / 2;
579 :
580 :
581 0 : for (Int kk = 0; kk < nBeamChans; ++kk) {
582 0 : Double fratio = beamFreqs(kk) / beamFreqs(nBeamChans - 1);
583 : //cerr << "fratio " << fratio << endl;
584 0 : Float convRatio = convSamp / convSampling;
585 0 : blcin[3] = kk;
586 0 : trcin[3] = kk;
587 : //Slicer slout(blcin, trcin, Slicer::endIsLast);
588 0 : Matrix<Complex> o1(i1.shape(), Complex(0.0));
589 0 : Matrix<Complex> o2(i2.shape(), Complex(0.0));
590 0 : for (Int yy = 0; yy < i1.shape()(1); ++yy) {
591 : //Int nyy= (Double(yy-cenY)*fratio) + cenY;
592 0 : Double nyy = (Double((yy - cenY) * convRatio) / fratio) + cenY;
593 0 : Double cyy = ceil(nyy);
594 0 : Double fyy = floor(nyy);
595 0 : Int iy = nyy > fyy + 0.5 ? Int(cyy) : Int(fyy);
596 0 : if (cyy < 2 * cenY && fyy >= 0.0)
597 0 : for (Int xx = 0; xx < i1.shape()(0); ++xx) {
598 : //Int nxx= Int(Double(xx-cenX)*fratio) + cenX;
599 0 : Double nxx = Int(Double((xx - cenX) * convRatio) / fratio) + cenX;
600 0 : Double cxx = ceil(nxx);
601 0 : Double fxx = floor(nxx);
602 0 : Int ix = nxx > fxx + 0.5 ? Int(cxx) : Int(fxx);
603 0 : if (cxx < 2 * cenX && fxx >= 0.0) {
604 : //Double dist=sqrt((nxx-cxx)*(nxx-cxx)+(nyy-cyy)*(nyy-cyy))/sqrt(2.0);
605 : //o1(xx, yy)=float(1-dist)*i1(fxx, fyy)+ dist*i1(cxx,cyy);
606 0 : o1(xx, yy) = i1(ix, iy);
607 : //o2(xx, yy)=i2(nxx, nyy);
608 : //o2(xx, yy)=float(1-dist)*i2(fxx, fyy)+ dist*i2(cxx,cyy);
609 0 : o2(xx, yy) = i2(ix, iy);
610 : }
611 : }
612 : }
613 0 : twoDPB.putSlice(o1.reform(elshape), blcin);
614 0 : twoDPB2.putSlice(o2.reform(elshape), blcin);
615 0 : }
616 :
617 0 : }
618 :
619 : /*
620 : {
621 : TempImage<Float> screen2(TiledShape(pbShape, IPosition(4, pbShape(0), pbShape(1), 1, 1)), coords, memtobeused);
622 : // Matrix<Float> screenoo(convSize_p, convSize_p);
623 : //screenoo.set(1.0);
624 : //screen2.putSlice(screenoo,start);
625 : //screen2.set(1.0);
626 : for (Int k=0; k < nBeamChans; ++k){
627 : blcin[3]=k;
628 : trcin[3]=k;
629 : Slicer slin(blcin, trcin, Slicer::endIsLast);
630 : SubImage<Float> subim(screen2, slin, true);
631 : SubImage<Complex> subout(twoDPB2, slin, true);
632 : subim.set(1.0);
633 : //twoDPB.putSlice(screen, start);
634 : sj_p->applySquare(subim, subim, vb, 0);
635 : //// LatticeExpr<Complex> le(subim);
636 : //// subout.copyData(le);
637 : ///// LatticeFFT::cfft2d(subout);
638 :
639 : LatticeFFT::rcfft(subout, subim, true, false);
640 : IPosition iblc(4, 0, 3*subout.shape()(1)/8, 0, 0);
641 : IPosition itrc(4, 0, 5*subout.shape()(1)/8, 0, 0);
642 : for(Int x=subout.shape()(0)/2; x <(5*subout.shape()(0)/8); ++x){
643 :
644 : iblc[0]=x-subout.shape()(0)/2;
645 : itrc[0]=x-subout.shape()(0)/2;
646 : Slicer isl(iblc, itrc, Slicer::endIsLast);
647 : iblc[0]=x;
648 : subout.putSlice(subout.getSlice(isl), iblc);
649 : }
650 : for(Int x=subout.shape()(0)/2+1; x <(5*subout.shape()(0)/8); ++x){
651 :
652 : iblc[0]=x;
653 : itrc[0]=x;
654 : Slicer isl(iblc, itrc, Slicer::endIsLast);
655 : iblc[0]=subout.shape()(0)-x;
656 : subout.putSlice(subout.getSlice(isl), iblc);
657 : if(x==(subout.shape()(0)-1)){
658 : iblc[0]=0;
659 : subout.putSlice(subout.getSlice(isl), iblc);
660 : }
661 : }
662 :
663 : }
664 :
665 : //sj_p->applySquare(screen2, screen2, vb, 0);
666 : //LatticeExpr<Complex> le(screen2);
667 : //twoDPB2.copyData(le);
668 : }
669 :
670 : */
671 :
672 : /*
673 : if(0) {
674 : CoordinateSystem ftCoords(coords);
675 : //directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
676 : //AlwaysAssert(directionIndex>=0, AipsError);
677 : dc=coords.directionCoordinate(directionIndex);
678 : Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
679 : Vector<Int> shape(2); shape(0)=convSize_p;shape(1)=convSize_p;
680 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
681 : //ftCoords.replaceCoordinate(*ftdc, directionIndex);
682 : delete ftdc; ftdc=0;
683 : ostringstream os1;
684 : os1 << "Screen_" << vb.fieldId() ;
685 : PagedImage<Complex> thisScreen(twoDPB2.shape(), ftCoords, String(os1));
686 : //LatticeExpr<Float> le(abs(twoDPB2));
687 : thisScreen.copyData(twoDPB2);
688 : LatticeFFT::cfft2d(thisScreen, false);
689 : PagedImage<Complex> thisScreen0(twoDPB.shape(), ftCoords, String("PB_")+String(os1));
690 : thisScreen0.copyData(twoDPB);
691 : LatticeFFT::cfft2d(thisScreen0, false);
692 : }
693 : */
694 : /*
695 : // Now FFT and get the result back
696 : //LatticeFFT::cfft2d(twoDPB);
697 : //LatticeFFT::cfft2d(twoDPB2);
698 :
699 : // Write out FT of screen as an image
700 : if(0) {
701 : CoordinateSystem ftCoords(coords);
702 : directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
703 : AlwaysAssert(directionIndex>=0, AipsError);
704 : dc=coords.directionCoordinate(directionIndex);
705 : Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
706 : Vector<Int> shape(2); shape(0)=convSize_p;shape(1)=convSize_p;
707 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
708 : ftCoords.replaceCoordinate(*ftdc, directionIndex);
709 : delete ftdc; ftdc=0;
710 : ostringstream os1;
711 : os1 << "FTScreen_" << vb.fieldId() ;
712 : PagedImage<Float> thisScreen(pbShape, ftCoords, String(os1));
713 : LatticeExpr<Float> le(abs(twoDPB2));
714 : thisScreen.copyData(le);
715 : }
716 : */
717 : //cerr << "twoDPB shape " << twoDPB.shape() << " slice shape " << IPosition(4, tempConvSize, tempConvSize, 1, 1) << endl;
718 0 : convFunc_p = twoDPB.getSlice(IPosition(4, 0, 0, 0, 0), IPosition(4, tempConvSize, tempConvSize, 1, 1), true);
719 0 : weightConvFunc_p.resize();
720 0 : weightConvFunc_p = twoDPB2.getSlice(IPosition(4, 0, 0, 0, nBeamChans - 1), IPosition(4, tempConvSize, tempConvSize, 1, 1), true);
721 : //convFunc/=max(abs(convFunc));
722 0 : Float maxAbsConvFunc = max(amplitude(weightConvFunc_p));
723 :
724 0 : Float minAbsConvFunc = min(amplitude(weightConvFunc_p));
725 : //cerr << "min max " << minAbsConvFunc << " " << maxAbsConvFunc << endl;
726 0 : convSupport_p = -1;
727 0 : Bool found = false;
728 : //Bool found2=true;
729 : //Int trial2=0;
730 0 : Int trial = 0;
731 0 : for (trial = tempConvSize / 2 - 2;trial > 0;trial--) {
732 : //Searching down a diagonal
733 0 : if (abs(weightConvFunc_p(tempConvSize / 2 - trial, tempConvSize / 2 - trial)) > (1.0e-3 * maxAbsConvFunc)) {
734 0 : found = true;
735 0 : trial = Int(sqrt(2.0 * Float(trial * trial)));
736 0 : break;
737 : }
738 : }
739 0 : if (!found) {
740 0 : if ((maxAbsConvFunc - minAbsConvFunc) > (1.0e-3 * maxAbsConvFunc))
741 0 : found = true;
742 : // if it drops by more than 2 magnitudes per pixel
743 0 : trial = (tempConvSize > (10 * convSampling)) ? 5 * convSampling : (tempConvSize / 2 - 4 * convSampling);
744 : }
745 :
746 0 : if (trial < 5 * convSampling)
747 0 : trial = (tempConvSize > (10 * convSampling)) ? 5 * convSampling : (tempConvSize / 2 - 4 * convSampling);
748 :
749 0 : if (found) {
750 0 : convSupport_p = Int(0.5 + Float(trial) / Float(convSampling)) + 1;
751 : }
752 : else {
753 : os << "Convolution function is misbehaved - support seems to be zero\n"
754 : << "Reasons can be: \n(1)The image definition not covering one or more of the pointings selected\n"
755 : << "(2) No unflagged data in a given pointing\n"
756 : << "(3) The entries in the POINTING subtable do not match the field being imaged."
757 : << "Please check, and try again with an empty POINTING subtable.)\n"
758 0 : << LogIO::EXCEPTION;
759 : }
760 :
761 : // Normalize such that plane 0 sums to 1 (when jumping in
762 : // steps of convSampling)
763 :
764 0 : Double pbSum = 0.0;
765 : //Double pbSum2=0.0;
766 :
767 :
768 :
769 0 : for (Int iy = -convSupport_p;iy <= convSupport_p;iy++) {
770 0 : for (Int ix = -convSupport_p;ix <= convSupport_p;ix++) {
771 0 : Complex val = convFunc_p(ix * convSampling + tempConvSize / 2,
772 0 : iy * convSampling + tempConvSize / 2);
773 0 : pbSum += real(val);
774 :
775 : //pbSum+=sqrt(real(val)*real(val)+ imag(val)*imag(val));
776 : }
777 : }
778 :
779 :
780 : //pbSum=sum(amplitude(convFunc_p))/Double(convSampling)/Double(convSampling);
781 :
782 0 : if (pbSum > 0.0) {
783 0 : convFunc_p *= Complex(1.0 / pbSum, 0.0);
784 : }
785 : else {
786 : os << "Convolution function integral is not positive"
787 0 : << LogIO::EXCEPTION;
788 : }
789 :
790 : //##########################################
791 : os << "Convolution support = " << convSupport_p
792 : << " pixels in Fourier plane"
793 0 : << LogIO::POST;
794 :
795 0 : convSupportBlock_p.resize(actualConvIndex_p + 1);
796 0 : convSizes_p.resize(actualConvIndex_p + 1);
797 : //Only one beam for now...but later this should be able to
798 : // take all the beams for the different antennas.
799 0 : convSupportBlock_p[actualConvIndex_p] = new Vector<Int>(1);
800 0 : convSizes_p[actualConvIndex_p] = new Vector<Int>(1);
801 0 : (*(convSupportBlock_p[actualConvIndex_p]))[0] = convSupport_p;
802 0 : convFunctions_p.resize(actualConvIndex_p + 1);
803 0 : convWeights_p.resize(actualConvIndex_p + 1);
804 0 : convFunctions_p[actualConvIndex_p] = new Array<Complex>();
805 0 : convWeights_p[actualConvIndex_p] = new Array<Complex>();
806 0 : Int newConvSize = 2 * (convSupport_p + 2) * convSampling;
807 : //NEED to chop this right ...and in the centre
808 0 : if (newConvSize >= tempConvSize)
809 0 : newConvSize = tempConvSize;
810 :
811 0 : IPosition blc(4, (tempConvSize / 2) - (newConvSize / 2),
812 0 : (tempConvSize / 2) - (newConvSize / 2), 0, 0);
813 0 : IPosition trc(4, (tempConvSize / 2) + (newConvSize / 2 - 1),
814 0 : (tempConvSize / 2) + (newConvSize / 2 - 1), 0, nBeamChans - 1);
815 0 : convFunctions_p[actualConvIndex_p]->resize(IPosition(5, newConvSize, newConvSize, 1, nBeamChans, 1));
816 : //cerr << "convFunc shape " << (convFunctions_p[actualConvIndex_p])->shape() <<
817 : //" " << " twoDPB shape " <<twoDPB.get(false)(blc,trc).shape() << endl;
818 0 : convFunctions_p[actualConvIndex_p]->copyMatchingPart(twoDPB.get(false)(blc, trc));//*Complex(1.0/pbSum,0.0));
819 0 : convSize_p = newConvSize;
820 0 : convWeights_p[actualConvIndex_p]->resize(IPosition(5, newConvSize, newConvSize, 1, nBeamChans, 1));
821 0 : convWeights_p[actualConvIndex_p]->copyMatchingPart(twoDPB2.get(false)(blc, trc));//*Complex(1.0/pbSum2,0.0));
822 0 : blc.resize(5, false);
823 0 : trc.resize(5, false);
824 0 : blc = IPosition(5, 0, 0, 0, 0, 0);
825 0 : trc = IPosition(5, newConvSize - 1, newConvSize - 1, 0, 0, 0);
826 0 : for (Int bc = 0; bc < nBeamChans; ++bc) {
827 0 : blc[3] = bc;
828 0 : trc[3] = bc;
829 0 : pbSum = real(sum((*convFunctions_p[actualConvIndex_p])(blc, trc))) / Double(convSampling) / Double(convSampling);
830 0 : (*convFunctions_p[actualConvIndex_p])(blc, trc) = (*convFunctions_p[actualConvIndex_p])(blc, trc) / pbSum;
831 0 : (*convWeights_p[actualConvIndex_p])(blc, trc) = (*convWeights_p[actualConvIndex_p])(blc, trc) / pbSum;
832 : }
833 :
834 :
835 :
836 0 : convFunc_p.resize();//break any reference
837 0 : (*convSizes_p[actualConvIndex_p])[0] = convSize_p;
838 0 : doneMainConv_p[actualConvIndex_p] = true;
839 :
840 0 : }
841 : else {
842 0 : convSize_p = (*convSizes_p[actualConvIndex_p])[0];
843 :
844 : }
845 :
846 : //Apply the shift phase gradient
847 0 : convFunc.resize();
848 0 : weightConvFunc.resize();
849 0 : convFunc.assign(*(convFunctions_p[actualConvIndex_p]));
850 0 : weightConvFunc.assign(*(convWeights_p[actualConvIndex_p]));
851 : Bool copyconv, copywgt;
852 0 : Complex* cv = convFunc.getStorage(copyconv);
853 0 : Complex* wcv = weightConvFunc.getStorage(copywgt);
854 : //cerr << "Field " << vb.fieldId() << " spw " << vb.spectralWindow() << " phase grad: " << pixFieldDir << endl;
855 :
856 0 : for (Int nc = 0; nc < nBeamChans; ++nc) {
857 0 : Int planeoffset = nc * convSize_p * convSize_p;
858 0 : for (Int iy = 0;iy < convSize_p;iy++) {
859 : Double cy, sy;
860 : Int offset;
861 :
862 0 : SINCOS(Double(iy - convSize_p / 2) * pixFieldDir(1), sy, cy);
863 0 : Complex phy(cy, sy);
864 0 : offset = iy * convSize_p + planeoffset;
865 0 : for (Int ix = 0;ix < convSize_p;ix++) {
866 : Double cx, sx;
867 0 : SINCOS(Double(ix - convSize_p / 2) * pixFieldDir(0), sx, cx);
868 0 : Complex phx(cx, sx);
869 0 : cv[ix + offset] = cv[ix + offset] * phx * phy;
870 0 : wcv[ix + offset] = wcv[ix + offset] * phx * phy;
871 :
872 : }
873 : }
874 : }
875 0 : convFunc.putStorage(cv, copyconv);
876 0 : weightConvFunc.putStorage(wcv, copywgt);
877 0 : convsize.resize();
878 0 : convsize = *(convSizes_p[actualConvIndex_p]);
879 0 : convSupport.resize();
880 0 : convSupport = (*(convSupportBlock_p[actualConvIndex_p]));
881 :
882 :
883 0 : }
884 :
885 0 : void SimplePBConvFunc::setSkyJones(SkyJones* sj) {
886 0 : sj_p = sj;
887 0 : }
888 :
889 0 : void SimplePBConvFunc::findUsefulChannels(std::vector<double>& freqs, const vi::VisBuffer2& vb) {
890 0 : Int spw = vb.spectralWindows()(0);
891 : //bandName_p=vb.subtableColumns().spectralWindow().name()(spw);
892 0 : Vector<Double> spwfreq = vb.subtableColumns().spectralWindow().chanFreq()(spw);
893 :
894 :
895 0 : double tol = (max(spwfreq)) * 1.0 / 100.0;
896 0 : Double spwfreqwidth = abs(Vector<Double>(vb.subtableColumns().spectralWindow().chanWidth()(spw))(0));
897 0 : if (tol < spwfreqwidth)
898 0 : tol = spwfreqwidth;
899 0 : Double topFreq = max(spwfreq);
900 0 : Double bottomFreq = min(spwfreq);
901 0 : uint nchan = std::round(topFreq - bottomFreq) / tol;
902 0 : if (nchan == 0)
903 0 : nchan = 1;
904 :
905 : ///////////////TESTOO
906 : //nchan = 1;
907 : //tol = (topFreq - bottomFreq) / 2.0;
908 : //////////////////////////
909 0 : freqs.resize(nchan);
910 0 : for (uint k = 0; k < nchan; ++k) {
911 0 : freqs[k] = k * tol + bottomFreq;
912 : }
913 :
914 0 : }
915 :
916 0 : void SimplePBConvFunc::findUsefulChannels(std::vector<double>& freqs, const vi::VisBuffer2& vb, const std::pair<double, double>& frange) {
917 :
918 0 : freqs.resize(0);
919 0 : std::vector<double> freqsforvb;
920 0 : findUsefulChannels(freqsforvb, vb);
921 : double fmin, fmax;
922 0 : if (std::get<0>(frange) > std::get<1>(frange)) {
923 0 : fmin = std::get<1>(frange);
924 0 : fmax = std::get<0>(frange);
925 : }
926 : else {
927 0 : fmin = std::get<0>(frange);
928 0 : fmax = std::get<1>(frange);
929 : }
930 0 : for (auto f : freqsforvb) {
931 0 : if (f >= fmin && f <= fmax)
932 0 : freqs.push_back(f);
933 : }
934 :
935 :
936 0 : }
937 :
938 :
939 :
940 :
941 :
942 0 : void SimplePBConvFunc::findUsefulChannels(Vector<Int>& chanMap, Vector<Double>& chanFreqs, const vi::VisBuffer2& vb, const Vector<Double>& freq) {
943 :
944 0 : Int spw = vb.spectralWindows()(0);
945 0 : bandName_p = vb.subtableColumns().spectralWindow().name()(spw);
946 0 : Vector<Double> spwfreq = vb.subtableColumns().spectralWindow().chanFreq()(spw);
947 0 : Double spwfreqwidth = abs(Vector<Double>(vb.subtableColumns().spectralWindow().chanWidth()(spw))(0));
948 0 : chanMap.resize(freq.nelements());
949 0 : Vector<Double> localfreq = vb.getFrequencies(0, MFrequency::TOPO);
950 0 : Double minfreq = min(freq);
951 0 : Double maxfreq = max(freq);
952 0 : Double origwidth = freq.nelements() == 1 ? 1e12 : (maxfreq - minfreq) / (freq.nelements() - 1);
953 : ///Fractional bandwidth which will trigger mutiple PB in one spw
954 :
955 0 : Double tol = (max(spwfreq)) * 0.5 / 100;
956 0 : if (tol < origwidth / 2.0) tol = origwidth / 2.0;
957 0 : Double topFreq = max(spwfreq);
958 0 : while (topFreq > maxfreq) {
959 0 : topFreq -= tol;
960 : }
961 : // For large tol
962 0 : if (topFreq < minfreq)
963 0 : topFreq += tol;
964 : //Int nchan=Int(lround((max(freq)-min(freq))/tol));
965 0 : Double bottomFreq = topFreq;
966 0 : Int nchan = 0;
967 : //cerr << std::setprecision(12) << "top " << bottomFreq << " minFreq " << minfreq << " maxFreq " << maxfreq << endl;
968 0 : while (bottomFreq > minfreq) {
969 0 : ++nchan;
970 0 : bottomFreq -= tol;
971 : }
972 0 : if (nchan > 1) {
973 0 : nchan -= 1;
974 0 : bottomFreq += tol;
975 : }
976 :
977 : //cerr << "TOLERA " << tol << " nchan " << nchan << " bot " << bottomFreq << " vb.nchan " << vb.nChannels() << endl;
978 : //Number of beams that matters are the ones in the data
979 0 : if (nchan > vb.nChannels()) {
980 0 : nchan = vb.nChannels();
981 0 : tol = spwfreqwidth;
982 0 : bottomFreq = min(localfreq);
983 : }
984 :
985 0 : chanFreqs.resize();
986 0 : if (nchan >= (Int)(freq.nelements() - 1)) { indgen(chanMap); chanFreqs = freq; return; }
987 0 : if ((nchan == 0) || (freq.nelements() == 1)) { chanFreqs = Vector<Double>(1, bottomFreq);chanMap.set(0); return; }
988 :
989 0 : chanFreqs.resize(nchan);
990 0 : for (Int k = 0; k < nchan; ++k) {
991 0 : chanFreqs[k] = k * tol + bottomFreq;
992 : }
993 :
994 :
995 0 : Int activechan = 0;
996 0 : chanMap.set(0);
997 0 : for (uInt k = 0; k < chanMap.nelements(); ++k) {
998 0 : Double mindiff = DBL_MAX;
999 0 : Int nearestchan = 0;
1000 0 : while ((activechan < nchan) && Double(abs(freq[k] - chanFreqs[activechan])) > Double(tol / Double(2.0))) {
1001 0 : if (mindiff > Double(abs(freq[k] - chanFreqs[activechan]))) {
1002 0 : mindiff = Double(abs(freq[k] - chanFreqs[activechan]));
1003 0 : nearestchan = activechan;
1004 : }
1005 :
1006 : // cerr << "k " << k << " atcivechan " << activechan << " comparison "
1007 : // << freq[k] << " " << chanFreqs[activechan] << endl;
1008 0 : ++activechan;
1009 : }
1010 0 : if (activechan != nchan) {
1011 0 : chanMap[k] = activechan;
1012 : }
1013 : //////////////////
1014 : else {
1015 : //cerr << std::setprecision(10) << "freq diffs " << freq[k]-chanFreqs << " TOL " << tol/2.0 << endl;
1016 0 : chanMap[k] = nearestchan;
1017 : }
1018 : ///////////////////////////
1019 0 : activechan = 0;
1020 : }
1021 : //cerr << "chanfreqs " << chanFreqs << endl;
1022 : //cerr << "USEFULchan " << chanMap << endl;
1023 0 : return;
1024 0 : }
1025 : /*
1026 :
1027 : void SimplePBConvFunc::findUsefulChannels(Vector<Int>& chanMap, Vector<Double>& chanFreqs, const vi::VisBuffer2& vb, const Vector<Double>& freq){
1028 :
1029 :
1030 : Int spw=vb.spectralWindows()(0);
1031 : bandName_p=vb.subtableColumns().spectralWindow().name()(spw);
1032 : chanMap.resize(freq.nelements());
1033 : Vector<Double> localfreq=vb.getFrequencies(0, MFrequency::TOPO);
1034 : Double minfreq=min(freq);
1035 :
1036 : Double origwidth=freq.nelements()==1 ? 1e12 : (max(freq)-min(freq))/(freq.nelements()-1);
1037 : ///Fractional bandwidth which will trigger mutiple PB in one spw
1038 : Double tol=(max(freq))*0.5/100;
1039 :
1040 : Int nchan=Int(lround((max(freq)-min(freq))/tol));
1041 :
1042 : //cerr << "TOLERA " << tol << " nchan " << nchan << " vb.nchan " << vb.nChannel() << endl;
1043 : //Number of beams that matters are the ones in the data
1044 : if(nchan > vb.nChannels())
1045 : nchan=vb.nChannels();
1046 :
1047 : if(tol < origwidth) tol=origwidth;
1048 : chanFreqs.resize();
1049 : if(nchan >= (Int)(freq.nelements()-1)) { indgen(chanMap); chanFreqs=freq; return;}
1050 : if((nchan==0) || (freq.nelements()==1)) { chanFreqs=Vector<Double>(1, freq[0]);chanMap.set(0); return;}
1051 :
1052 : //readjust the tolerance...
1053 : tol=(max(freq)-min(freq)+origwidth)/Double(nchan);
1054 : chanFreqs.resize(nchan);
1055 : for (Int k=0; k < nchan; ++k)
1056 : chanFreqs[k]=minfreq-origwidth+tol/2.0+tol*Double(k);
1057 : Int activechan=0;
1058 : chanMap.set(0);
1059 : for (uInt k=0; k < chanMap.nelements(); ++k){
1060 : Double mindiff=DBL_MAX;
1061 : Int nearestchan=0;
1062 : while((activechan< nchan) && Double(abs(freq[k]-chanFreqs[activechan])) > Double(tol/Double(2.0))){
1063 : if(mindiff > Double(abs(freq[k]-chanFreqs[activechan]))){
1064 : mindiff=Double(abs(freq[k]-chanFreqs[activechan]));
1065 : nearestchan=activechan;
1066 : }
1067 :
1068 : // cerr << "k " << k << " atcivechan " << activechan << " comparison "
1069 : // << freq[k] << " " << chanFreqs[activechan] << endl;
1070 : ++activechan;
1071 : }
1072 : if(activechan != nchan){
1073 : chanMap[k]=activechan;
1074 : }
1075 : //////////////////
1076 : else{
1077 : //cerr << std::setprecision(10) << "freq diffs " << freq[k]-chanFreqs << " TOL " << tol/2.0 << endl;
1078 : chanMap[k]=nearestchan;
1079 : }
1080 : ///////////////////////////
1081 : activechan=0;
1082 : }
1083 : //cerr << chanMap << endl;
1084 : return;
1085 : }
1086 : */
1087 0 : Bool SimplePBConvFunc::checkPBOfField(const vi::VisBuffer2& vb) {
1088 : //Int fieldid=vb.fieldId();
1089 0 : String msid = vb.msName(true);
1090 : /*
1091 : if(convFunctionMap_p.ndefined() > 0){
1092 : if (((fluxScale_p.shape()[3] != nchan_p) || (fluxScale_p.shape()[2] != npol_p)) && calcFluxScale_p){
1093 : convFunctionMap_p.clear();
1094 : }
1095 : }
1096 : // if you rename the ms might be a problem
1097 : String mapid=msid+String("_")+String::toString(fieldid);
1098 : if(convFunctionMap_p.ndefined() == 0){
1099 : convFunctionMap_p.define(mapid, 0);
1100 : actualConvIndex_p=0;
1101 : if(calcFluxScale_p){
1102 : // 0ne channel only is needed to keep track of pb coverage
1103 : if(fluxScale_p.shape().nelements()==0){
1104 : fluxScale_p=TempImage<Float>(IPosition(4,nx_p,ny_p,npol_p,1), csys_p);
1105 : fluxScale_p.set(0.0);
1106 : }
1107 : filledFluxScale_p=false;
1108 : }
1109 : return false;
1110 : }
1111 :
1112 : if(!convFunctionMap_p.isDefined(mapid)){
1113 : actualConvIndex_p=convFunctionMap_p.ndefined();
1114 : convFunctionMap_p.define(mapid, actualConvIndex_p);
1115 : return false;
1116 : }
1117 : else{
1118 : actualConvIndex_p=convFunctionMap_p(mapid);
1119 : convFunc_p.resize(); // break any reference
1120 : weightConvFunc_p.resize();
1121 : //Here we will need to use the right xyPlane for different PA range.
1122 : //convFunc_p.reference(convFunctions_p[actualConvIndex_p]->xyPlane(0));
1123 : //weightConvFunc_p.reference(convWeights_p[actualConvIndex_p]->xyPlane(0));
1124 : //Again this for one time of antenna only later should be fixed for all
1125 : // antennas independently
1126 : convSupport_p=(*convSupportBlock_p[actualConvIndex_p])[0];
1127 : convSize_p=(*convSizes_p[actualConvIndex_p])[0];
1128 :
1129 : }
1130 : */
1131 :
1132 0 : return true;
1133 :
1134 :
1135 :
1136 0 : }
1137 :
1138 :
1139 0 : ImageInterface<Float>& SimplePBConvFunc::getFluxScaleImage() {
1140 :
1141 0 : if (!calcFluxScale_p)
1142 0 : throw(AipsError("Programmer error: Cannot get flux scale"));
1143 0 : if (!filledFluxScale_p) {
1144 0 : IPosition blc = fluxScale_p.shape();
1145 0 : IPosition trc = fluxScale_p.shape();
1146 0 : blc(0) = 0; blc(1) = 0; trc(0) = nx_p - 1; trc(1) = ny_p - 1;
1147 :
1148 0 : for (Int j = 0; j < fluxScale_p.shape()(2); ++j) {
1149 0 : for (Int k = 0; k < fluxScale_p.shape()(3); ++k) {
1150 :
1151 0 : blc(2) = j; trc(2) = j;
1152 0 : blc(3) = k; trc(3) = k;
1153 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1154 0 : SubImage<Float> fscalesub(fluxScale_p, sl, true);
1155 : Float planeMax;
1156 0 : LatticeExprNode LEN = max(fscalesub);
1157 0 : planeMax = LEN.getFloat();
1158 0 : if (planeMax != 0) {
1159 0 : fscalesub.copyData((LatticeExpr<Float>) (fscalesub / planeMax));
1160 :
1161 : }
1162 0 : }
1163 : }
1164 : /*
1165 : if(0) {
1166 : ostringstream os2;
1167 : os2 << "ALL_" << "BEAMS" ;
1168 : PagedImage<Float> thisScreen2(fluxScale_p.shape(), fluxScale_p.coordinates(), String(os2));
1169 : thisScreen2.copyData(fluxScale_p);
1170 : }
1171 : */
1172 :
1173 0 : filledFluxScale_p = true;
1174 0 : }
1175 :
1176 :
1177 0 : return fluxScale_p;
1178 :
1179 : }
1180 :
1181 :
1182 0 : Bool SimplePBConvFunc::toRecord(RecordInterface& rec) {
1183 0 : Int numConv = convFunctions_p.nelements();
1184 : // not saving the protected variables as they are generated by
1185 : // the first call to storeImageParams
1186 : try {
1187 0 : rec.define("name", "SimplePBConvFunc");
1188 0 : rec.define("numconv", numConv);
1189 : //cerr << "num of conv " << numConv << " " << convFunctionMap_p.ndefined() << " " <<convFunctions_p.nelements() << endl;
1190 0 : std::map<String, Int>::iterator it = vbConvIndex_p.begin();
1191 0 : for (Int k = 0; k < numConv; ++k) {
1192 0 : rec.define("convfunctions" + String::toString(k), *(convFunctions_p[k]));
1193 0 : rec.define("convweights" + String::toString(k), *(convWeights_p[k]));
1194 0 : rec.define("convsizes" + String::toString(k), *(convSizes_p[k]));
1195 0 : rec.define("convsupportblock" + String::toString(k), *(convSupportBlock_p[k]));
1196 : //cerr << "k " << k << " key " << convFunctionMap_p.getKey(k) << " val " << convFunctionMap_p.getVal(k) << endl;
1197 0 : rec.define(String("key") + String::toString(k), it->first);
1198 0 : rec.define(String("val") + String::toString(k), it->second);
1199 0 : it++;
1200 : }
1201 0 : rec.define("pbclass", Int(pbClass_p));
1202 0 : rec.define("actualconvindex", actualConvIndex_p);
1203 0 : rec.define("donemainconv", doneMainConv_p);
1204 0 : rec.define("usepointingtable", usePointingTable_p);
1205 : //The following is not needed ..can be regenerated
1206 : //rec.define("pointingpix", pointingPix_p);
1207 : }
1208 0 : catch (AipsError& x) {
1209 0 : return false;
1210 0 : }
1211 0 : return true;
1212 : }
1213 :
1214 0 : Bool SimplePBConvFunc::fromRecord(String& err, const RecordInterface& rec, Bool calcFluxneeded) {
1215 0 : Int numConv = 0;
1216 : //make sure storeImageParams is triggered
1217 0 : nchan_p = 0;
1218 :
1219 : try {
1220 0 : if (!rec.isDefined("name") || rec.asString("name") != "SimplePBConvFunc") {
1221 0 : throw(AipsError("Wrong record to recover SimplePBConvFunc from"));
1222 : }
1223 0 : rec.get("numconv", numConv);
1224 0 : convFunctions_p.resize(numConv, true, false);
1225 0 : convSupportBlock_p.resize(numConv, true, false);
1226 0 : convWeights_p.resize(numConv, true, false);
1227 0 : convSizes_p.resize(numConv, true, false);
1228 0 : convFunctionMap_p.clear();
1229 0 : vbConvIndex_p.erase(vbConvIndex_p.begin(), vbConvIndex_p.end());
1230 0 : for (Int k = 0; k < numConv; ++k) {
1231 0 : convFunctions_p[k] = new Array<Complex>();
1232 0 : convWeights_p[k] = new Array<Complex>();
1233 0 : convSizes_p[k] = new Vector<Int>();
1234 0 : convSupportBlock_p[k] = new Vector<Int>();
1235 0 : rec.get("convfunctions" + String::toString(k), *(convFunctions_p[k]));
1236 0 : rec.get("convsupportblock" + String::toString(k), *(convSupportBlock_p[k]));
1237 0 : rec.get("convweights" + String::toString(k), *(convWeights_p[k]));
1238 0 : rec.get("convsizes" + String::toString(k), *(convSizes_p[k]));
1239 0 : String key;
1240 : Int val;
1241 0 : rec.get(String("key") + String::toString(k), key);
1242 0 : rec.get(String("val") + String::toString(k), val);
1243 0 : vbConvIndex_p[key] = val;
1244 0 : ant1PointVal_p.clear();
1245 0 : ant1PointingCache_p.resize();
1246 : //convFunctionMap_p.define(key,val);
1247 0 : }
1248 0 : pbClass_p = static_cast<PBMathInterface::PBClass>(rec.asInt("pbclass"));
1249 0 : rec.get("actualconvindex", actualConvIndex_p);
1250 0 : rec.get("usepointingtable", usePointingTable_p);
1251 0 : pointingPix_p.resize();
1252 : //rec.get("pointingpix", pointingPix_p);
1253 0 : calcFluxScale_p = calcFluxneeded;
1254 :
1255 : }
1256 0 : catch (AipsError& x) {
1257 0 : err = x.getMesg();
1258 0 : return false;
1259 0 : }
1260 0 : return true;
1261 :
1262 : }
1263 0 : void SimplePBConvFunc::addPBToFlux(const vi::VisBuffer2& vb) {
1264 0 : if (calcFluxScale_p) {
1265 0 : if (fluxScale_p.shape().nelements() == 0)
1266 : { //cerr << "nx_p " << nx_p << endl;
1267 0 : calcFluxScale_p = False;
1268 0 : return;
1269 : }
1270 0 : Vector<Int> pixdepoint(2, -100000);
1271 0 : convertArray(pixdepoint, thePix_p);
1272 0 : if ((pixdepoint(0) >= 0) && (pixdepoint(0) < pointingPix_p.shape()[0]) && (pixdepoint(1) >= 0) && (pixdepoint(1) < pointingPix_p.shape()[1]) && !pointingPix_p(pixdepoint(0), pixdepoint(1))) {
1273 0 : TempImage<Float> thispb(fluxScale_p.shape(), fluxScale_p.coordinates());
1274 0 : thispb.set(1.0);
1275 0 : sj_p->applySquare(thispb, thispb, vb, 0);
1276 0 : LatticeExpr<Float> le(fluxScale_p + thispb);
1277 0 : fluxScale_p.copyData(le);
1278 0 : pointingPix_p(pixdepoint(0), pixdepoint(1)) = true;
1279 : //LatticeExprNode LEN = max(fluxScale_p);
1280 : //Float maxsca=LEN.getFloat();
1281 : //Tempporary fix when cubesky is chunking...do not add on
1282 : //already defined position
1283 : //if(maxsca > 1.98){
1284 : // cerr << "avoiding subtract " << endl;
1285 : //fluxScale_p.copyData(LatticeExpr<Float>(fluxScale_p-thispb));
1286 :
1287 : //}
1288 : /*
1289 : if(0) {
1290 : ostringstream os1;
1291 : os1 << "SINGLE_" << vb.fieldId() ;
1292 : PagedImage<Float> thisScreen(fluxScale_p.shape(), fluxScale_p.coordinates(), String(os1));
1293 : thisScreen.copyData(thispb);
1294 : ostringstream os2;
1295 : os2 << "ALL_" << vb.fieldId() ;
1296 : PagedImage<Float> thisScreen2(fluxScale_p.shape(), fluxScale_p.coordinates(), String(os2));
1297 : thisScreen2.copyData(fluxScale_p);
1298 : }
1299 : */
1300 0 : }
1301 0 : }
1302 :
1303 : }
1304 :
1305 0 : void SimplePBConvFunc::sliceFluxScale(Int npol) {
1306 0 : IPosition fshp = fluxScale_p.shape();
1307 0 : if (fshp(2) > npol) {
1308 0 : npol_p = npol;
1309 : // use first npol planes...
1310 0 : IPosition blc(4, 0, 0, 0, 0);
1311 0 : IPosition trc(4, fluxScale_p.shape()(0) - 1, fluxScale_p.shape()(1) - 1, npol - 1, fluxScale_p.shape()(3) - 1);
1312 0 : Slicer sl = Slicer(blc, trc, Slicer::endIsLast);
1313 : //writeable if possible
1314 0 : SubImage<Float> fluxScaleSub = SubImage<Float>(fluxScale_p, sl, true);
1315 0 : fluxScale_p = TempImage<Float>(fluxScaleSub.shape(), fluxScaleSub.coordinates());
1316 0 : LatticeExpr<Float> le(fluxScaleSub);
1317 0 : fluxScale_p.copyData(le);
1318 0 : }
1319 0 : }
1320 0 : void SimplePBConvFunc::setAWConvFuncHolder(std::shared_ptr<AWConvFuncHolder> awptr){
1321 0 : awConvs_p=awptr;
1322 :
1323 0 : }
1324 0 : std::shared_ptr<AWConvFuncHolder> SimplePBConvFunc::getAWConvFuncHolder(){
1325 0 : return awConvs_p;
1326 : }
1327 :
1328 : } //# END of name space REFIM
1329 : } //# NAMESPACE CASA - END
1330 :
1331 :
1332 :
1333 :
1334 :
1335 :
1336 :
|
