Line data Source code
1 : /*
2 : * MSUVBin.cc implementation of gridding MSs to a gridded MS
3 : //#
4 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
5 :
6 : //# Copyright (C) 2014-2015
7 : //# Associated Universities, Inc. Washington DC, USA.
8 : //#
9 : //# This library is free software; you can redistribute it and/or modify it
10 : //# under the terms of the GNU General Public License as published by
11 : //# the Free Software Foundation; either version 3 of the License, or (at your
12 : //# option) any later version.
13 : //#
14 : //# This library is distributed in the hope that it will be useful, but WITHOUT
15 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
17 : //# License for more details.
18 : //#
19 : //# You should have received a copy of the GNU General Public License
20 : //# along with this library; if not, write to the Free Software Foundation,
21 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
22 : //#
23 :
24 : //# Postal address:
25 : //# National Radio Astronomy Observatory
26 : //# 520 Edgemont Road
27 : //# Charlottesville, VA 22903-2475 USA
28 : //#
29 : *
30 : * Created on: Feb 4, 2014
31 : * Author: kgolap
32 : */
33 :
34 : #include <casacore/casa/Arrays/ArrayMath.h>
35 : #include <casacore/casa/Arrays/Array.h>
36 : #include <casacore/casa/Arrays/Matrix.h>
37 : #include <casacore/casa/Arrays/Cube.h>
38 : #include <casacore/casa/Arrays/Vector.h>
39 : #include <casacore/casa/OS/HostInfo.h>
40 :
41 : #include <casacore/casa/System/ProgressMeter.h>
42 : #include <casacore/casa/Quanta/QuantumHolder.h>
43 : #include <casacore/casa/Utilities/CompositeNumber.h>
44 : #include <casacore/measures/Measures/MeasTable.h>
45 : #include <casacore/ms/MeasurementSets/MSPolColumns.h>
46 : #include <casacore/ms/MeasurementSets/MSPolarization.h>
47 : #include <mstransform/MSTransform/MSUVBin.h>
48 : #include <mstransform/MSTransform/MSTransformDataHandler.h>
49 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
50 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
51 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
52 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
53 : #include <casacore/images/Images/PagedImage.h>
54 : #include <casacore/images/Images/TempImage.h>
55 : #include <msvis/MSVis/MSUtil.h>
56 : #include <msvis/MSVis/VisBuffer.h>
57 : #include <msvis/MSVis/VisBuffer2Adapter.h>
58 : #include <imageanalysis/Utilities/SpectralImageUtil.h>
59 : #include <msvis/MSVis/VisibilityIterator2.h>
60 : #include <casacore/scimath/Mathematics/FFTPack.h>
61 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
62 : #include <wcslib/wcsconfig.h> /** HAVE_SINCOS **/
63 : #include <math.h>
64 : #ifdef _OPENMP
65 : #include <omp.h>
66 : #endif
67 : #if HAVE_SINCOS
68 : #define SINCOS(a,s,c) sincos(a,&s,&c)
69 : #else
70 : #define SINCOS(a,s,c) \
71 : s = sin(a); \
72 : c = cos(a)
73 : #endif
74 :
75 :
76 : typedef unsigned long long ooLong;
77 :
78 : using namespace casacore;
79 : namespace casa { //# NAMESPACE CASA - BEGIN
80 :
81 0 : MSUVBin::MSUVBin():nx_p(0), ny_p(0), nchan_p(0), npol_p(0),existOut_p(false),numVis_p(), sumWeight_p(){
82 0 : outMsPtr_p=NULL;
83 0 : outMSName_p="OutMS.ms";
84 0 : memFraction_p=0.5;
85 0 : }
86 0 : MSUVBin::MSUVBin(const MDirection& phaseCenter,
87 0 : const Int nx, const Int ny, const Int nchan, const Int npol, Quantity cellx, Quantity celly, Quantity freqStart, Quantity freqStep, Float memFraction, Bool dow, Bool doflag):
88 0 : existOut_p(false)
89 :
90 : {
91 0 : phaseCenter_p=phaseCenter;
92 0 : nx_p=nx;
93 0 : ny_p=ny;
94 0 : nchan_p=nchan;
95 0 : npol_p=npol;
96 0 : numVis_p.resize(npol, nchan);
97 0 : sumWeight_p.resize(npol, nchan);
98 0 : numVis_p.set(0);
99 0 : sumWeight_p.set(0);
100 0 : deltas_p.resize(2);
101 0 : deltas_p[0]=cellx.getValue("rad");
102 0 : deltas_p[1]=celly.getValue("rad");
103 0 : freqStart_p=freqStart.getValue("Hz");
104 0 : freqStep_p=freqStep.getValue("Hz");
105 0 : outMsPtr_p=NULL;
106 0 : outMSName_p="OutMS.ms";
107 0 : memFraction_p=memFraction;
108 0 : doW_p=dow;
109 0 : doFlag_p=doflag;
110 :
111 0 : }
112 :
113 0 : MSUVBin::~MSUVBin(){
114 : //close the measurementsets
115 0 : for (uInt k=0; k < mss_p.nelements(); ++k){
116 0 : (const_cast<MeasurementSet *>(mss_p[k]))->unlock();
117 0 : *(const_cast<MeasurementSet *>(mss_p[k]))=MeasurementSet();
118 : }
119 0 : }
120 0 : Bool MSUVBin::selectData(const String& msname, const String& spw, const String& field,
121 : const String& baseline, const String& scan,
122 : const String& uvrange, const String& taql,
123 : const String& subarray, const String& correlation,const String& intent, const String& obs){
124 :
125 0 : Vector<Int> fakestep = Vector<Int> (1, 1);
126 :
127 0 : MSTransformDataHandler mshandler(msname, doFlag_p ? Table::Update: Table::Old);
128 : //No very well documented what the step is supposed to do
129 : //using the default value here
130 0 : if(mshandler.setmsselect(spw, field, baseline, scan, uvrange,
131 : taql, fakestep, subarray, correlation, intent, obs)){
132 0 : mss_p.resize(mss_p.nelements()+1);
133 0 : mshandler.makeSelection();
134 : // have to make a copy here as the selected ms pointer of mshandler
135 : // dies with mshandler
136 0 : mss_p[mss_p.nelements()-1]=new MeasurementSet(*mshandler.getSelectedInputMS());
137 : }
138 : else
139 0 : return false;
140 :
141 : // cerr << "num of ms" << mss_p.nelements() << " ms0 " << (mss_p[0])->tableName() << " nrows " << (mss_p[0])->nrow() << endl;
142 0 : return true;
143 0 : }
144 :
145 0 : void MSUVBin::setOutputMS(const String& msname){
146 0 : outMSName_p=msname;
147 0 : checkOutputGridParams();
148 0 : }
149 :
150 :
151 0 : void MSUVBin::createOutputMS(const Int nrrow){
152 0 : if(Table::isReadable(outMSName_p)){
153 0 : Int oldnrows=recoverGridInfo(outMSName_p);
154 0 : if(oldnrows != nrrow)
155 0 : throw(AipsError("Number of grid points requested "+ String::toString(nrrow)+ " does not match "+String::toString(oldnrows)+ " in outputms"));
156 0 : existOut_p=true;
157 0 : return;
158 : }
159 0 : if(mss_p.nelements()==0)
160 0 : throw(AipsError("no ms selected for input yet"));
161 0 : Vector<Int> tileShape(3);
162 0 : tileShape[0]=4; tileShape[1]=200; tileShape[2]=500;
163 0 : auto wtSpec = MSColumns(*mss_p[0]).weightSpectrum();
164 : //auto createdWeightSpectrumCols = !wtSpec.isNull() && wtSpec.isDefined(0);
165 0 : outMsPtr_p = MSTransformDataHandler::setupMS(outMSName_p, nchan_p, npol_p,
166 0 : Vector<MS::PredefinedColumns>(1, MS::DATA),
167 0 : True, tileShape);
168 0 : outMsPtr_p->addRow(nrrow, true);
169 : //cerr << "mss Info " << mss_p[0]->tableName() << " " << mss_p[0]->nrow() <<endl;
170 0 : MSTransformDataHandler::addOptionalColumns(mss_p[0]->spectralWindow(),
171 0 : outMsPtr_p->spectralWindow());
172 : ///Setup pointing why this is not done in setupMS
173 : {
174 0 : SetupNewTable pointingSetup(outMsPtr_p->pointingTableName(),
175 0 : MSPointing::requiredTableDesc(), Table::New);
176 : // POINTING can be large, set some sensible defaults for storageMgrs
177 0 : IncrementalStMan ismPointing ("ISMPointing");
178 0 : StandardStMan ssmPointing("SSMPointing", 32768);
179 0 : pointingSetup.bindAll(ismPointing, true);
180 0 : pointingSetup.bindColumn(MSPointing::columnName(MSPointing::DIRECTION),
181 : ssmPointing);
182 0 : pointingSetup.bindColumn(MSPointing::columnName(MSPointing::TARGET),
183 : ssmPointing);
184 0 : pointingSetup.bindColumn(MSPointing::columnName(MSPointing::TIME),
185 : ssmPointing);
186 0 : outMsPtr_p->rwKeywordSet().defineTable(MS::keywordName(MS::POINTING),
187 0 : Table(pointingSetup));
188 0 : outMsPtr_p->initRefs();
189 0 : }
190 0 : TableCopy::copySubTables(outMsPtr_p->pointing(), (mss_p[0])->pointing());
191 0 : if(doW_p){
192 0 : TiledShapeStMan tiledStMan("TiledImagWgtSpectrum", tileShape);
193 0 : outMsPtr_p->addColumn(ArrayColumnDesc<Float>("IMAGINARY_WEIGHT_SPECTRUM",2), tiledStMan);
194 0 : ArrayColumn<Float> imagWgtSpecCol(*outMsPtr_p, "IMAGINARY_WEIGHT_SPECTRUM");
195 0 : imagWgtSpecCol.fillColumn(Matrix<Float>(npol_p, nchan_p, 0.0));
196 0 : }
197 0 : MSColumns msc(*outMsPtr_p);
198 0 : msc.data().fillColumn(Matrix<Complex>(npol_p, nchan_p, Complex(0.0)));
199 0 : msc.flagRow().fillColumn(true);
200 0 : msc.flag().fillColumn(Matrix<Bool>(npol_p, nchan_p, true));
201 0 : msc.weight().fillColumn(Vector<Float>(npol_p, 0.0));
202 0 : msc.sigma().fillColumn(Vector<Float>(npol_p, 0.0));
203 0 : msc.weightSpectrum().fillColumn(Matrix<Float>(npol_p, nchan_p, 0.0));
204 0 : msc.antenna1().fillColumn(0);
205 0 : msc.antenna2().fillColumn(0);
206 : //change array id for every 2000 rows
207 0 : Vector<Int> arrayID(nrrow,0);
208 0 : Int blk=2000;
209 0 : Int naid=nrrow/blk;
210 0 : IPosition blc(1,0);
211 0 : IPosition trc(1,0);
212 0 : for (Int k=0; k < naid; ++k){
213 0 : blc[0]=k*blk;
214 0 : trc[0]= (k < (naid-1)) ? (k+1)*blk-1 : (k+1)*blk+nrrow%blk-1;
215 0 : arrayID(blc, trc)=k;
216 : }
217 : //cerr << "MINMAX array ID " << min(arrayID) << " " << max(arrayID) << endl;
218 0 : msc.arrayId().putColumn(arrayID);
219 : // Note: we suspect this flush with sync=true (fsync!) can be removed (it is
220 : // commented out in other functions in this same file. See CAS-11946.
221 0 : outMsPtr_p->flush(true);
222 :
223 0 : existOut_p=false;
224 :
225 0 : }
226 0 : bool MSUVBin::checkOutputGridParams(){
227 : //Table does not exist nothing to check
228 0 : if(!Table::isReadable(outMSName_p))
229 0 : return true;
230 0 : MeasurementSet theGridMS(outMSName_p, Table::Old);
231 0 : if(!theGridMS.keywordSet().isDefined("MSUVBIN")){
232 0 : throw(AipsError("The ms "+outMSName_p+" was not made with the UV binner"));
233 : }
234 0 : Record rec=theGridMS.keywordSet().asRecord("MSUVBIN");
235 : Int nx, ny, nchan, npol;
236 0 : rec.get("nx", nx);
237 0 : rec.get("ny", ny);
238 0 : rec.get("nchan", nchan);
239 0 : rec.get("npol", npol);
240 0 : if(!doFlag_p){
241 0 : if(nx != nx_p || ny != ny_p || nchan != nchan_p || npol != npol_p){
242 0 : LogIO os(LogOrigin("MSUVBin", "checkGridParams", WHERE));
243 0 : os << LogIO::SEVERE << "Output is a binned ms of "<< nx << " by " << ny << " with "<< npol << " pol and "<< nchan << " channels where as user has given [" << nx_p << ", " << ny_p << ", " << npol_p << ", " << nchan_p << "]" << LogIO::POST;
244 0 : return false;
245 0 : }
246 : }
247 : CoordinateSystem *tempCoordsys;
248 0 : tempCoordsys=CoordinateSystem::restore(rec, "csys");
249 0 : if(tempCoordsys==NULL)
250 0 : throw(AipsError("could recover grid info from ms"));
251 :
252 0 : Vector<Double> pixelPhaseCenter(2);
253 0 : pixelPhaseCenter(0) = Double( nx / 2 );
254 0 : pixelPhaseCenter(1) = Double( ny / 2 );
255 0 : MDirection phcen;
256 0 : (tempCoordsys->directionCoordinate(0)).toWorld(phcen, pixelPhaseCenter);
257 0 : if(!doFlag_p && ((phaseCenter_p.getRefPtr()->getType()) != ((phcen.getRefPtr())->getType()))){
258 0 : throw(AipsError("Requested output frame is not the same as what was used with existant "+outMSName_p));
259 : }
260 0 : MVDirection mvphcen=phcen.getValue();
261 0 : Double sep=mvphcen.separation(phaseCenter_p.getValue());
262 0 : if(!doFlag_p && (sep > 2.0*(tempCoordsys->directionCoordinate(0)).increment()(0)))
263 : {
264 0 : std::ostringstream oss;
265 0 : oss << "stored is " << phcen.toString();
266 0 : oss << " user requested " <<phaseCenter_p.toString();
267 0 : throw(AipsError("phasecentre requested is not the same as stored in "+outMSName_p+ " " +String(oss.str())));
268 0 : }
269 :
270 : // When transferring flags we don't care about the input params we take
271 : // what is stored
272 0 : if(doFlag_p){
273 0 : nx_p=nx;
274 0 : ny_p=ny;
275 0 : nchan_p=nchan;
276 0 : npol_p=npol;
277 0 : csys_p=*tempCoordsys;
278 0 : phaseCenter_p=phcen;
279 : }
280 :
281 0 : delete tempCoordsys;
282 0 : return true;
283 :
284 0 : }
285 0 : Int MSUVBin::recoverGridInfo(const String& msname){
286 0 : outMsPtr_p=new MeasurementSet(msname, Table::Update);
287 0 : if(!outMsPtr_p->keywordSet().isDefined("MSUVBIN")){
288 0 : throw(AipsError("The ms "+msname+" was not made with the UV binner"));
289 : }
290 0 : Record rec=outMsPtr_p->keywordSet().asRecord("MSUVBIN");
291 0 : rec.get("nx", nx_p);
292 0 : rec.get("ny", ny_p);
293 0 : rec.get("nchan", nchan_p);
294 0 : rec.get("npol", npol_p);
295 0 : LogIO os(LogOrigin("MSUVBin", "recoverInfo", WHERE));
296 0 : os << "Output is a binned ms of "<< nx_p << " by " << ny_p << " with "<< npol_p << " pol and "<< nchan_p << " channel " << endl;
297 : CoordinateSystem *tempCoordsys;
298 0 : tempCoordsys=CoordinateSystem::restore(rec, "csys");
299 0 : if(tempCoordsys==NULL)
300 0 : throw(AipsError("could recover grid info from ms"));
301 0 : csys_p=*tempCoordsys;
302 0 : delete tempCoordsys;
303 0 : numVis_p.resize();
304 0 : rec.get("numvis", numVis_p);
305 0 : sumWeight_p.resize();
306 0 : rec.get("sumweight", sumWeight_p);
307 0 : return outMsPtr_p->nrow();
308 :
309 :
310 0 : }
311 :
312 0 : void MSUVBin::storeGridInfo(){
313 : //cerr << "numVis " << numVis_p << endl;
314 : //cerr << "sumWgt " << sumWeight_p << endl;
315 0 : if(existOut_p){
316 0 : outMsPtr_p->rwKeywordSet().asrwRecord("MSUVBIN").define("numvis", numVis_p);
317 0 : outMsPtr_p->rwKeywordSet().asrwRecord("MSUVBIN").define("sumweight", sumWeight_p);
318 0 : return;
319 : }
320 0 : Record rec;
321 0 : rec.define("nx", nx_p);
322 0 : rec.define("ny", ny_p);
323 0 : rec.define("nchan", nchan_p);
324 0 : rec.define("npol", npol_p);
325 0 : rec.define("numvis", numVis_p);
326 0 : rec.define("sumweight", sumWeight_p);
327 0 : csys_p.save(rec, "csys");
328 0 : outMsPtr_p->rwKeywordSet().defineRecord("MSUVBIN", rec);
329 :
330 :
331 0 : }
332 0 : void MSUVBin::setTileCache()
333 : {
334 0 : if(outMsPtr_p->tableType() ==Table::Memory){
335 0 : return;
336 : }
337 0 : const ColumnDescSet & cds = outMsPtr_p->tableDesc ().columnDescSet ();
338 :
339 : //uInt startrow = 0;
340 :
341 0 : Vector<String> columns (6);
342 0 : columns (0) = MS::columnName (MS::DATA); // complex
343 0 : columns (1) = MS::columnName (MS::FLAG); // boolean
344 0 : columns (2) = MS::columnName (MS::WEIGHT_SPECTRUM); // float
345 0 : columns (3) = MS::columnName (MS::WEIGHT); // float
346 0 : columns (4) = MS::columnName (MS::SIGMA); // float
347 0 : columns (5) = MS::columnName (MS::UVW); // double
348 :
349 0 : for (uInt k = 0; k < columns.nelements (); ++k) {
350 :
351 0 : if (! cds.isDefined (columns (k)))
352 : {
353 0 : continue;
354 : }
355 :
356 : try {
357 : //////////////////
358 : //////Temporary fix for virtual ms of multiple real ms's ...miracle of statics
359 : //////setting the cache size of hypercube at row 0 of each ms.
360 : ///will not work if each subms of a virtual ms has multi hypecube being
361 : ///accessed.
362 :
363 : {
364 0 : Bool usesTiles=false;
365 0 : String dataManType = RODataManAccessor (*outMsPtr_p, columns[k], true).dataManagerType ();
366 : //cerr << "column " << columns[k] << " dataman "<< dataManType << endl;
367 0 : usesTiles = dataManType.contains ("Tiled");
368 0 : if(usesTiles){
369 0 : ROTiledStManAccessor tacc (*outMsPtr_p, columns[k], true);
370 0 : tacc.clearCaches (); //One tile only for now ...seems to work faster
371 :
372 :
373 :
374 :
375 :
376 : /// If some bucketSize is 0...there is trouble in setting cache
377 : /// but if slicer is used it gushes anyways if one does not set cache
378 : /// need to fix the 0 bucket size in the filler anyways...then this is not needed
379 :
380 :
381 0 : tacc.setCacheSize (0, 1);
382 :
383 :
384 0 : }
385 0 : }
386 : }
387 0 : catch (AipsError x) {
388 : // cerr << "Data man type " << dataManType << " " << dataManType.contains ("Tiled") << " && " << (!String (cdesc.dataManagerGroup ()).empty ()) << endl;
389 : // cerr << "Failed to set settilecache due to " << x.getMesg () << " column " << columns[k] <<endl;
390 : //It failed so leave the caching as is
391 0 : continue;
392 0 : }
393 : }
394 :
395 0 : }
396 :
397 0 : Bool MSUVBin::fillOutputMS(){
398 : // Double imagevol=Double(nx_p)*Double(ny_p)*Double(npol_p)*Double(nchan_p)*12.0/1e6; //in MB
399 : // Double memoryMB=Double(HostInfo::memoryFree())/1024.0 ;
400 : Bool retval;
401 0 : retval= fillNewBigOutputMS();
402 0 : return retval;
403 : }
404 :
405 0 : Bool MSUVBin::fillNewBigOutputMS(){
406 0 : if(mss_p.nelements()==0)
407 0 : throw(AipsError("No valid input MSs defined"));
408 0 : Vector<Double> incr;
409 0 : Vector<Int> cent;
410 0 : Matrix<Double> uvw;
411 0 : Double numOfFlagsBefore=0;
412 0 : Double numOfFlagsAfter=0;
413 : //need to build or recover csys at this stage
414 0 : makeCoordsys();
415 0 : Double reffreq=SpectralImageUtil::worldFreq(csys_p, Double(nchan_p/2));
416 0 : Int nrrows=makeUVW(reffreq, incr, cent, uvw);
417 0 : Vector<Bool> rowFlag(nrrows);
418 0 : Vector<Int> ant1(nrrows);
419 0 : Vector<Int> ant2(nrrows);
420 0 : Vector<Double> timeCen(nrrows);
421 0 : Matrix<Float> wght(npol_p, nrrows);
422 0 : createOutputMS(nrrows);
423 0 : setTileCache();
424 0 : MSColumns msc(*outMsPtr_p);
425 0 : vi::VisibilityIterator2 iter(mss_p, vi::SortColumns(), doFlag_p);
426 0 : vi::VisBuffer2* vb=iter.getVisBuffer();
427 0 : iter.originChunks();
428 0 : iter.origin();
429 0 : if(existOut_p){
430 0 : msc.weight().getColumn(wght);
431 0 : msc.flagRow().getColumn(rowFlag);
432 0 : msc.antenna1().getColumn(ant1);
433 0 : msc.antenna2().getColumn(ant2);
434 0 : msc.time().getColumn(timeCen);
435 0 : msc.uvw().getColumn(uvw);
436 :
437 : }
438 : else{
439 0 : wght.set(0.0);
440 0 : rowFlag.set(true);
441 0 : timeCen.set(vb->time()(0));
442 : }
443 : ////////////////////////////////////////////////
444 0 : Cube<Complex> convFunc;
445 0 : Vector<Int> convSupport;
446 : Double wScale;
447 : Int convSampling, convSize;
448 0 : if(doW_p){
449 0 : makeWConv(iter, convFunc, convSupport, wScale, convSampling, convSize);
450 : //makeSFConv(convFunc, convSupport, wScale, convSampling, convSize);
451 : // cerr << "convSupport0 " << convFunc.shape() << " " << convSupport << endl;
452 : }
453 : /////////////////////////////////////////////////
454 0 : Int usableNchan=Int(Double(HostInfo::memoryFree())*memFraction_p*1024.0/Double(npol_p)/Double(nx_p*ny_p)/(doW_p ? 16.0: 12.0));
455 0 : if(usableNchan < nchan_p)
456 0 : cerr << "Maximum nchan per pass " << min(usableNchan, nchan_p) << endl;
457 0 : Int npass=nchan_p%usableNchan==0 ? nchan_p/usableNchan : nchan_p/usableNchan+1;
458 0 : if(npass >1)
459 0 : cerr << "Due to lack of memory will be doing " << npass << " passes through data"<< endl;
460 0 : for (Int pass=0; pass < npass; ++pass){
461 0 : Int startchan=pass*usableNchan;
462 0 : Int endchan=pass==nchan_p/usableNchan ? nchan_p%usableNchan+ startchan-1 : (pass+1)*usableNchan -1 ;
463 :
464 0 : Cube<Complex> grid(npol_p, endchan-startchan+1, nrrows);
465 : //cerr << "shape " << grid.shape() << endl;
466 0 : Cube<Complex> wghtSpec;
467 0 : Cube<Float> realWghtSpec(npol_p, endchan-startchan+1, nrrows);
468 0 : Cube<Float> imagWghtSpec;
469 0 : Cube<Bool> flag(npol_p, endchan-startchan+1, nrrows);
470 : //Matrix<Int> locuv;
471 :
472 :
473 0 : iter.originChunks();
474 0 : iter.origin();
475 0 : vbutil_p=VisBufferUtil(*vb);
476 0 : if(existOut_p){
477 : //recover the previous data for summing
478 0 : Slicer elslice(IPosition(2, 0, startchan), IPosition(2,npol_p, endchan-startchan+1));
479 0 : msc.data().getColumn(elslice, grid);
480 0 : msc.weightSpectrum().getColumn(elslice, realWghtSpec);
481 0 : msc.flag().getColumn(elslice, flag);
482 0 : wghtSpec.resize(realWghtSpec.shape());
483 0 : setReal(wghtSpec, realWghtSpec);
484 0 : if(outMsPtr_p->tableDesc().columnDescSet().isDefined("IMAGINARY_WEIGHT_SPECTRUM")){
485 0 : ArrayColumn<Float> imagWgtSpecCol(*outMsPtr_p, "IMAGINARY_WEIGHT_SPECTRUM");
486 0 : imagWghtSpec.resize(realWghtSpec.shape());
487 0 : imagWgtSpecCol.getColumn(elslice, imagWghtSpec);
488 0 : setImag(wghtSpec, imagWghtSpec);
489 0 : }
490 0 : if(doW_p){
491 0 : imagWghtSpec.resize();
492 0 : realWghtSpec.resize();
493 : }
494 : //multiply the data with weight here
495 : {
496 0 : for (Int iz=0; iz< grid.shape()(2); ++iz){
497 0 : for(Int iy=0; iy < grid.shape()(1); ++iy){
498 0 : for(Int ix=0; ix < grid.shape()(0); ++ix){
499 0 : grid(ix,iy,iz)= grid(ix,iy,iz)*wghtSpec(ix,iy,iz);
500 : }
501 : }
502 : }
503 : }
504 0 : if(!doW_p)
505 0 : wghtSpec.resize();
506 :
507 :
508 0 : }
509 : else{
510 0 : grid.set(Complex(0));
511 0 : if(doW_p){
512 0 : wghtSpec.resize(realWghtSpec.shape());
513 0 : wghtSpec.set(0);
514 0 : realWghtSpec.resize();
515 : }
516 0 : realWghtSpec.set(0);
517 0 : flag.set(true);
518 :
519 : //cerr << "Zeroing grid " << grid.shape() << endl;
520 : //outMsPtr_p->addRow(nrrows, true);
521 : }
522 0 : if(npass > 1)
523 0 : cerr <<"Pass " << pass ;
524 : ProgressMeter pm(1.0, Double(nrrows),"Gridding data",
525 0 : "", "", "", true);
526 0 : Double rowsDone=0.0;
527 : //cerr << "Before: Num of flagged model " << ntrue(uvw.row(2) ==(-666.0)) << endl;
528 :
529 0 : for (iter.originChunks(); iter.moreChunks(); iter.nextChunk()){
530 0 : for(iter.origin(); iter.more(); iter.next()){
531 0 : if(doFlag_p){
532 0 : numOfFlagsBefore += ntrue(vb->flagCube());
533 : //cerr << " before " << ntrue(vb->flagCube()) << endl;
534 0 : Cube<Bool> datFlag=vb->flagCube();
535 0 : locateFlagFromGrid(*vb, datFlag,
536 : realWghtSpec,
537 : flag, rowFlag, uvw, ant1,
538 : ant2, timeCen, startchan, endchan);
539 : //cerr << " after " << ntrue(datFlag) << endl;
540 0 : numOfFlagsAfter += ntrue(datFlag);
541 0 : iter.writeFlag(datFlag);
542 :
543 0 : }
544 : else{
545 0 : if(doW_p){
546 : //gridDataConv(*vb, grid, wght, wghtSpec,flag, rowFlag, uvw,
547 : // ant1,ant2,timeCen, startchan, endchan, convFunc, convSupport, wScale, convSampling);
548 0 : gridDataConvThr(*vb, grid, wghtSpec,flag, rowFlag, uvw,
549 : ant1,ant2,timeCen, startchan, endchan, convFunc, convSupport, wScale, convSampling);
550 :
551 : }
552 : else{
553 0 : gridData(*vb, grid, wght, realWghtSpec,flag, rowFlag, uvw,
554 : ant1,ant2,timeCen, startchan, endchan);
555 : }
556 :
557 :
558 : }
559 0 : rowsDone+=Double(vb->nRows());
560 :
561 :
562 :
563 0 : pm.update(rowsDone);
564 :
565 :
566 : }
567 :
568 : }
569 :
570 0 : imagWghtSpec.resize();
571 0 : if(doW_p){
572 0 : realWghtSpec.resize(wghtSpec.shape());
573 0 : realWghtSpec=real(wghtSpec);
574 0 : imagWghtSpec.resize(wghtSpec.shape());
575 0 : imagWghtSpec=imag(wghtSpec);
576 : }
577 : else{
578 0 : wghtSpec.resize(realWghtSpec.shape());
579 0 : wghtSpec.set(0.0);
580 0 : setReal(wghtSpec, realWghtSpec);
581 : }
582 :
583 : //Weight Correct the data
584 : {
585 0 : for (Int iz=0; iz< grid.shape()(2); ++iz){
586 : /* Float medweight= 0.0;
587 : if(max(abs(wghtSpec.xyPlane(iz))) > 0.0){
588 : medweight=median(wghtSpec.xyPlane(iz));
589 : medweight
590 : }
591 : */
592 : /* if(max(wghtSpec.xyPlane(iz)) > 0.0)
593 : cerr << iz << " median " << median(wghtSpec.xyPlane(iz)) << " min " << min(wghtSpec.xyPlane(iz)) << " max " << max(wghtSpec.xyPlane(iz)) << endl;
594 : */
595 :
596 0 : for(Int iy=0; iy < grid.shape()(1); ++iy){
597 0 : for(Int ix=0; ix < grid.shape()(0); ++ix){
598 0 : if(!flag(ix,iy,iz) && wghtSpec(ix,iy, iz) > 0.0){
599 0 : grid(ix,iy, iz)=grid(ix,iy,iz)/wghtSpec(ix,iy,iz);
600 : }
601 : else{
602 0 : grid(ix,iy,iz)=0.0;
603 0 : wghtSpec(ix,iy,iz)=0.0;
604 : }
605 :
606 : //grid(ix,iy,iz)= (!flag(ix,iy,iz) && fabs(wghtSpec(ix,iy, iz)) > 1e-3*medweight) ? grid(ix,iy,iz)/wghtSpec(ix,iy,iz): Complex(0);
607 : }
608 : }
609 :
610 :
611 : } // iz
612 : }
613 :
614 : //cerr << "After: Num of flagged model " << ntrue(uvw.row(2) ==(-666.0)) << endl;
615 0 : saveData(grid, flag, rowFlag, realWghtSpec, uvw, ant1, ant2, timeCen, startchan, endchan, imagWghtSpec);
616 0 : }
617 0 : if(doFlag_p){
618 0 : LogIO os(LogOrigin("MSUVBin", "TransferFlags", WHERE));
619 0 : os << "Number of flags before transfer " << numOfFlagsBefore << "\n";
620 0 : os << "Number of flags after transfer " << numOfFlagsAfter << LogIO::POST;
621 0 : }
622 0 : storeGridInfo();
623 0 : return true;
624 0 : }
625 :
626 0 : Bool MSUVBin::fillSmallOutputMS(){
627 0 : if(mss_p.nelements()==0)
628 0 : throw(AipsError("No valid input MSs defined"));
629 0 : Vector<Double> incr;
630 0 : Vector<Int> cent;
631 0 : Matrix<Double> uvw;
632 : //need to build or recover csys at this stage
633 0 : makeCoordsys();
634 0 : Double reffreq=SpectralImageUtil::worldFreq(csys_p, Double(nchan_p/2));
635 0 : Int nrrows=makeUVW(reffreq, incr, cent, uvw);
636 0 : Cube<Complex> grid(npol_p, nchan_p, nrrows);
637 0 : Matrix<Float> wght(npol_p, nrrows);
638 0 : Cube<Float> wghtSpec(npol_p, nchan_p, nrrows);
639 0 : Cube<Bool> flag(npol_p, nchan_p, nrrows);
640 0 : Vector<Bool> rowFlag(nrrows);
641 0 : Vector<Int> ant1(nrrows);
642 0 : Vector<Int> ant2(nrrows);
643 0 : Vector<Double> timeCen(nrrows);
644 0 : createOutputMS(nrrows);
645 0 : Matrix<Int> locuv;
646 0 : vi::VisibilityIterator2 iter(mss_p, vi::SortColumns(), false);
647 0 : vi::VisBuffer2* vb=iter.getVisBuffer();
648 :
649 0 : iter.originChunks();
650 0 : iter.origin();
651 0 : vbutil_p=VisBufferUtil(*vb);
652 0 : if(existOut_p){
653 : //recover the previous data for summing
654 :
655 0 : MSColumns msc(*outMsPtr_p);
656 0 : msc.data().getColumn(grid);
657 0 : msc.weightSpectrum().getColumn(wghtSpec);
658 0 : msc.weight().getColumn(wght);
659 0 : msc.flag().getColumn(flag);
660 0 : msc.flagRow().getColumn(rowFlag);
661 0 : msc.antenna1().getColumn(ant1);
662 0 : msc.antenna2().getColumn(ant2);
663 0 : msc.time().getColumn(timeCen);
664 0 : msc.uvw().getColumn(uvw);
665 0 : }
666 : else{
667 0 : grid.set(Complex(0));
668 0 : wght.set(0);
669 0 : wghtSpec.set(0);
670 0 : flag.set(true);
671 0 : rowFlag.set(true);
672 : //cerr << "SETTING time to val" << vb->time()(0) << endl;
673 0 : timeCen.set(vb->time()(0));
674 : //outMsPtr_p->addRow(nrrows, true);
675 : }
676 : ProgressMeter pm(1.0, Double(nrrows),"Gridding data",
677 0 : "", "", "", true);
678 0 : Double rowsDone=0.0;
679 0 : for (iter.originChunks(); iter.moreChunks(); iter.nextChunk()){
680 0 : for(iter.origin(); iter.more(); iter.next()){
681 0 : locateuvw(locuv, incr, cent, vb->uvw());
682 0 : gridData(*vb, grid, wght, wghtSpec,flag, rowFlag, uvw,
683 : ant1,ant2,timeCen, locuv);
684 0 : rowsDone+=Double(vb->nRows());
685 0 : pm.update(rowsDone);
686 : }
687 :
688 : }
689 :
690 0 : saveData(grid, flag, rowFlag, wghtSpec, wght, uvw, ant1, ant2, timeCen);
691 0 : storeGridInfo();
692 0 : return true;
693 0 : }
694 :
695 0 : Bool MSUVBin::fillBigOutputMS(){
696 :
697 0 : if(mss_p.nelements()==0)
698 0 : throw(AipsError("No valid input MSs defined"));
699 0 : Vector<Double> incr;
700 0 : Vector<Int> cent;
701 0 : Matrix<Double> uvw;
702 : //need to build or recover csys at this stage
703 0 : makeCoordsys();
704 0 : Double reffreq=SpectralImageUtil::worldFreq(csys_p, Double(nchan_p/2));
705 0 : Int nrrows=makeUVW(reffreq, incr, cent, uvw);
706 0 : createOutputMS(nrrows);
707 0 : if(!existOut_p)
708 0 : MSColumns(*outMsPtr_p).uvw().putColumn(uvw);
709 0 : nrrows=0;
710 0 : for (uInt k=0; k < mss_p.nelements() ; ++k)
711 0 : nrrows+=(mss_p[k])->nrow();
712 :
713 0 : vi::VisibilityIterator2 iter(mss_p, vi::SortColumns(), false);
714 0 : vi::VisBuffer2* vb=iter.getVisBuffer();
715 0 : iter.originChunks();
716 0 : iter.origin();
717 0 : vbutil_p=VisBufferUtil(*vb);
718 0 : Matrix<Int> locuv;
719 :
720 : ProgressMeter pm(1.0, Double(nrrows),"Gridding data",
721 0 : "", "", "", true);
722 0 : Double rowsDone=0.0;
723 0 : for (iter.originChunks(); iter.moreChunks(); iter.nextChunk()){
724 0 : for(iter.origin(); iter.more(); iter.next()){
725 :
726 : //locateuvw(locuv, incr, cent, vb->uvw());
727 0 : inplaceGridData(*vb);
728 : //
729 : //gridData(*vb, grid, wght, wghtSpec,flag, rowFlag, uvw,
730 : // ant1,ant2,timeCen, locuv);
731 0 : rowsDone+=Double(vb->nRows());
732 0 : pm.update(rowsDone);
733 : }
734 : }
735 :
736 0 : if(!existOut_p){
737 0 : fillSubTables();
738 : }
739 : ////Need to do the weight calculation here from spectral weight.
740 0 : weightSync();
741 0 : storeGridInfo();
742 0 : return true;
743 0 : }
744 :
745 0 : Int MSUVBin::makeUVW(const Double reffreq, Vector<Double>& increment,
746 : Vector<Int>& center, Matrix<Double>& uvw){
747 0 : Vector<Int> shp(2);
748 0 : shp(0)=nx_p; shp(1)=ny_p;
749 0 : Int directionIndex=csys_p.findCoordinate(Coordinate::DIRECTION);
750 0 : DirectionCoordinate thedir=csys_p.directionCoordinate(directionIndex);
751 0 : Coordinate *ftcoord=thedir.makeFourierCoordinate(Vector<Bool>(2, true), shp);
752 0 : increment=ftcoord->increment();
753 0 : increment *= C::c/reffreq;
754 : //Vector<Float> scale(2);
755 : //scale(0)=1.0/(nx_p*thedir.increment()(0));
756 : //scale(1)=1.0/(ny_p*thedir.increment()(1));
757 0 : center.resize(2);
758 0 : center(0)=nx_p/2;
759 0 : center(1)=ny_p/2;
760 0 : Vector<Int> npix(2);
761 0 : npix(0)=nx_p;
762 0 : npix(1)=ny_p;
763 0 : uvw.resize(3, nx_p*ny_p);
764 0 : uvw.set(-666.0);
765 0 : Vector<Double> px(2);
766 0 : Vector<Double> wld(2);
767 0 : Double lambd=C::c/reffreq;
768 0 : Int counter=0;
769 0 : for (Int k=0; k < ny_p;++k ){
770 0 : px(1)=Double(k);
771 : //px(1)=Double(k-ny_p/2)*C::c/reffreq*scale(1);
772 0 : for(Int j=0; j < nx_p; ++j){
773 : //px(0)=Double(j-nx_p/2)*C::c/reffreq*scale(0);
774 : //uvw(0, k*nx_p+j)=px(0);
775 : //uvw(1, k*nx_p+j)=px(1);
776 : //wld(0)=Double(j);
777 : //wld(1)=Double(k);
778 0 : px(0)=Double(j);
779 0 : if(ftcoord->toWorld(wld, px)){
780 : //cerr << "wld " << wld*C::c/reffreq << " px " << uvw(0,k*nx_p+j) << ", "<< uvw(1, k*nx_p+j) << endl;
781 0 : uvw(0, k*nx_p+j)=wld(0)*lambd;
782 0 : uvw(1, k*nx_p+j)=wld(1)*lambd;
783 0 : ++counter;
784 : }
785 : }
786 : }
787 0 : return counter;
788 0 : }
789 0 : void MSUVBin::makeCoordsys(){
790 :
791 : //if outMS already has the coordsys..it is recovered
792 0 : if(existOut_p)
793 0 : return;
794 : //
795 0 : Matrix<Double> xform(2,2);
796 0 : xform=0.0;xform.diagonal()=1.0;
797 0 : MVAngle ra=phaseCenter_p.getAngle("rad").getValue()(0);
798 0 : ra(0.0);
799 0 : MVAngle dec=phaseCenter_p.getAngle("rad").getValue()(1);
800 :
801 0 : DirectionCoordinate myRaDec(MDirection::Types(phaseCenter_p.getRefPtr()->getType()),
802 : Projection::SIN,
803 : ra.radian(), dec.radian(),
804 0 : deltas_p(0), deltas_p(1),
805 : xform,
806 0 : Double(nx_p)/2.0, Double(ny_p)/2.0);
807 :
808 : //Now to the spectral part
809 : //LSRK frame unless REST
810 0 : MFrequency::Types outFreqFrame=MFrequency::LSRK;
811 0 : Vector<MFrequency::Types> typesInData;
812 0 : MSUtil::getSpectralFrames(typesInData, *(mss_p[0]));
813 0 : if(anyEQ(typesInData, MFrequency::REST))
814 0 : outFreqFrame=MFrequency::REST;
815 : //make freqstart the centre of channel
816 0 : freqStart_p=freqStart_p+freqStep_p/2.0;
817 0 : SpectralCoordinate mySpectral(outFreqFrame, freqStart_p, freqStep_p, 0.0);
818 0 : MSColumns msc(*mss_p[0]);
819 0 : Int ddId=msc.dataDescId()(0);
820 0 : Int firstPolId=msc.dataDescription().polarizationId()(ddId);
821 0 : Int polType = Vector<Int>(msc.polarization().corrType()(firstPolId))(0);
822 0 : Bool isCircular=true;
823 0 : if (polType>=Stokes::XX && polType<=Stokes::YY)
824 0 : isCircular=false;
825 :
826 0 : if(isCircular){
827 0 : if(npol_p==4){
828 0 : whichStokes_p.resize(4);
829 :
830 0 : whichStokes_p(0)=Stokes::RR; whichStokes_p(1)=Stokes::RL;
831 0 : whichStokes_p(2)=Stokes::LR; whichStokes_p(3)=Stokes::LL;
832 : }
833 0 : else if(npol_p==2){
834 0 : whichStokes_p.resize(2);
835 0 : whichStokes_p(0)=Stokes::RR; whichStokes_p(1)=Stokes::LL;
836 : }
837 0 : else if(npol_p==1){
838 0 : whichStokes_p.resize(1);
839 0 : whichStokes_p(0)=Stokes::RR;
840 : }
841 : }
842 : else{
843 0 : if(npol_p==4){
844 0 : whichStokes_p.resize(4);
845 0 : whichStokes_p(0)=Stokes::XX; whichStokes_p(1)=Stokes::XY;
846 0 : whichStokes_p(2)=Stokes::YX; whichStokes_p(3)=Stokes::YY;
847 : }
848 0 : else if(npol_p==2){
849 0 : whichStokes_p.resize(2);
850 0 : whichStokes_p(0)=Stokes::XX; whichStokes_p(1)=Stokes::YY;
851 : }
852 0 : else if(npol_p==1){
853 0 : whichStokes_p.resize(1);
854 0 : whichStokes_p(0)=Stokes::XX;
855 : }
856 : }
857 0 : StokesCoordinate myStokes(whichStokes_p);
858 0 : csys_p=CoordinateSystem();
859 0 : csys_p.addCoordinate(myRaDec);
860 0 : csys_p.addCoordinate(myStokes);
861 0 : csys_p.addCoordinate(mySpectral);
862 0 : ObsInfo obsinf;
863 0 : obsinf.setObsDate(msc.timeMeas()(0));
864 0 : obsinf.setTelescope(msc.observation().telescopeName()(0));
865 0 : obsinf.setPointingCenter(phaseCenter_p.getValue());
866 0 : csys_p.setObsInfo(obsinf);
867 0 : }
868 :
869 0 : void MSUVBin::locateuvw(Matrix<Int>& locuv, const Vector<Double>& increment,
870 : const Vector<Int>& center, const Matrix<Double>& uvw){
871 0 : locuv.resize(2, uvw.shape()(1));
872 0 : for (Int k=0; k <uvw.shape()(1); ++k){
873 0 : for(Int j=0; j < 2; ++j)
874 0 : locuv(j,k)=Int(Double(center(j))+uvw(j,k)/increment(j)+0.5);
875 : }
876 :
877 :
878 0 : }
879 0 : Bool MSUVBin::datadescMap(const vi::VisBuffer2& vb, Double& fracbw){
880 0 : polMap_p.resize(vb.nCorrelations());
881 0 : polMap_p.set(-1);
882 : //Should ultimately find the real polarization in case it is a non common MS
883 0 : if(npol_p==4){
884 0 : if( polMap_p.nelements()==4)
885 0 : indgen(polMap_p);
886 0 : if(polMap_p.nelements()==2){
887 0 : polMap_p(0)=0;
888 0 : polMap_p(1)=3;
889 : }
890 : //This case sounds bad
891 0 : if(polMap_p.nelements()==1)
892 0 : polMap_p(0)=0;
893 : }
894 0 : if(npol_p==2){
895 0 : if( polMap_p.nelements()==4){
896 0 : polMap_p(0)=0;
897 0 : polMap_p(3)=1;
898 : }
899 0 : if(polMap_p.nelements()==2)
900 0 : indgen(polMap_p);
901 : //again a bad case ...not resolving this here
902 0 : if(polMap_p.nelements()==1)
903 0 : polMap_p(0)=0;
904 : }
905 0 : else if(npol_p==1){
906 0 : polMap_p(0)=0;
907 0 : if(polMap_p.nelements()==2)
908 0 : polMap_p(1)=0;
909 0 : if(polMap_p.nelements()==4)
910 0 : polMap_p(3)=0;
911 : }
912 : ////Now to chanmap
913 0 : fracbw=0.0;
914 0 : Vector<Double> f(1);
915 0 : Vector<Double> c(1,0.0);
916 0 : SpectralCoordinate spec=csys_p.spectralCoordinate(2);
917 0 : Vector<Double> visFreq=vb.getFrequencies(0, MFrequency::LSRK);
918 : //Vector<Double> visFreq;
919 : //vbutil_p.convertFrequency(visFreq, vb, MFrequency::LSRK);
920 0 : chanMap_p.resize(visFreq.nelements());
921 0 : chanMapRev_p.resize(nchan_p);
922 0 : for (Int chan=0; chan < nchan_p; ++chan){
923 0 : chanMapRev_p[chan].resize();
924 : }
925 0 : chanMap_p.set(-1);
926 0 : for (Int chan=0; chan < vb.nChannels(); ++chan){
927 0 : f[0]=visFreq[chan];
928 :
929 0 : if(spec.toPixel(c,f)){
930 :
931 0 : Int pixel=Int(floor(c(0)+0.5)); // round to chan freq at chan center
932 : //cout << " pixel " << pixel << " " << c << " f " << f << " chan " << chan << endl;
933 0 : if(pixel > -1 && pixel< nchan_p){
934 0 : chanMap_p(chan)=pixel;
935 0 : chanMapRev_p[pixel].resize(chanMapRev_p[pixel].nelements()+1, true);
936 0 : chanMapRev_p[pixel][chanMapRev_p[pixel].nelements()-1]=chan;
937 0 : c[0]=pixel;
938 0 : spec.toWorld(f, c);
939 : //cout << " pixel " << pixel << " chan " << chan << endl;
940 0 : if((abs(visFreq[chan]-f[0])/f[0]) > fracbw)
941 0 : fracbw=(abs(visFreq[chan]-f[0])/f[0]);
942 : }
943 : }
944 : }
945 :
946 : ///////////////////////
947 : //cout << "spw" << vb.spectralWindows()(0) << " rowid " << vb.rowIds()(0) << " max " << max(chanMap_p) << " sum "<< ntrue(chanMap_p > -1) << endl;
948 : /*if( ntrue(chanMap_p > -1) > 0){
949 : for (uInt k=0; k < visFreq.nelements() ; ++k){
950 : if(chanMap_p(k) > -1){
951 : cerr << " " <<k<< " "<< visFreq[k] << " = " << chanMap_p[k];
952 : }
953 : }
954 : cerr <<endl;
955 : cerr << "*************************************" << endl;
956 :
957 : }
958 : */
959 : ////////////////////////////////////
960 : //cerr << "spw " << vb.spectralWindows() << " fracbw " << fracbw << endl;
961 0 : if(allLT(chanMap_p ,0) || allLT(polMap_p , 0))
962 0 : return false;
963 :
964 0 : return true;
965 0 : }
966 :
967 0 : void MSUVBin::weightSync(){
968 0 : MSColumns msc(*outMsPtr_p);
969 0 : Matrix<Float> wght(npol_p, outMsPtr_p->nrow());
970 0 : Matrix<Float> wghtSpec;
971 0 : for (uInt row=0; row < outMsPtr_p->nrow(); ++row){
972 0 : for (uInt pol=0; pol < wght.shape()(0); ++pol){
973 :
974 0 : wghtSpec=msc.weightSpectrum().get(row);
975 : //cerr << "shape min max "<< median(wghtSpec.xyPlane(newrow).row(pol)) << " " << min(wghtSpec.xyPlane(newrow).row(pol)) << " "<< max(wghtSpec.xyPlane(newrow).row(pol)) << endl;
976 0 : wght(pol,row)=median(wghtSpec.row(pol));
977 : //cerr << "pol " << pol << " newrow "<< newrow << " weight "<< wght(pol, newrow) << endl;
978 : }
979 : }
980 0 : msc.weight().putColumn(wght);
981 0 : Matrix<Float> sigma(wght.shape());
982 0 : for (uInt k=0; k<wght.shape()[1]; ++k)
983 0 : for (uInt j=0; j<wght.shape()[0]; ++j)
984 0 : sigma(j,k)=wght(j,k)>0.0 ? 1/sqrt(wght(j,k)) :0.0;
985 0 : msc.sigma().putColumn(sigma);
986 :
987 :
988 0 : }
989 :
990 0 : void MSUVBin::inplaceGridData(const vi::VisBuffer2& vb){
991 : //we need polmap and chanmap;
992 : // no point of proceeding if nothing maps in this vb
993 : Double fracbw;
994 0 : if(!datadescMap(vb, fracbw)) return;
995 : //cerr << "fracbw " << fracbw << endl;
996 :
997 0 : if(fracbw >0.05)
998 0 : inplaceLargeBW(vb);
999 : else
1000 0 : inplaceSmallBW(vb);
1001 :
1002 :
1003 :
1004 :
1005 : }
1006 :
1007 0 : void MSUVBin::inplaceLargeBW(const vi::VisBuffer2& vb){
1008 : //Dang i thought the new vb will return Data or FloatData if correctedData was
1009 : //not there
1010 0 : Bool hasCorrected=!(MSMainColumns(vb.ms()).correctedData().isNull());
1011 0 : Int nrows=vb.nRows();
1012 : // Cube<Complex> grid(npol_p, nchan_p, nrows);
1013 : // Matrix<Float> wght(npol_p,nrows);
1014 : // Cube<Float> wghtSpec(npol_p,nchan_p,nrows);
1015 : // Cube<Bool> flag(npol_p, nchan_p,nrows);
1016 : // Vector<Bool> rowFlag(nrows);
1017 : // Matrix<Double> uvw(3, nrows);
1018 : // Vector<Int> ant1(nrows);
1019 : // Vector<Int> ant2(nrows);
1020 : // Vector<Double> timeCen(nrows);
1021 : // Vector<uInt> rowids(nrows);
1022 : //cerr << outMsPtr_p.null() << endl;
1023 0 : MSColumns msc(*outMsPtr_p);
1024 0 : SpectralCoordinate spec=csys_p.spectralCoordinate(2);
1025 0 : DirectionCoordinate thedir=csys_p.directionCoordinate(0);
1026 0 : Vector<Float> scale(2);
1027 0 : scale(0)=fabs(nx_p*thedir.increment()(0))/C::c;
1028 0 : scale(1)=fabs(ny_p*thedir.increment()(1))/C::c;
1029 : ///Index
1030 0 : Vector<Int> rowToIndex(nx_p*ny_p, -1);
1031 0 : Matrix<Int> locu(vb.nChannels(), nrows, -1);
1032 0 : Matrix<Int> locv(vb.nChannels(), nrows, -1);
1033 0 : Vector<Double> visFreq=vb.getFrequencies(0, MFrequency::LSRK);
1034 0 : Vector<Double> phasor;
1035 0 : Matrix<Double> eluvw;
1036 0 : Bool needRot=vbutil_p.rotateUVW(vb, phaseCenter_p, eluvw, phasor);
1037 0 : Int numUniq=0;
1038 0 : for (Int k=0; k <nrows; ++k){
1039 0 : for(Int chan=0; chan < vb.nChannels(); ++chan ){
1040 0 : if(chanMap_p(chan) >=0){
1041 0 : locv(chan, k)=Int(Double(ny_p)/2.0+eluvw(1,k)*visFreq(chan)*scale(1)+0.5);
1042 0 : locu(chan, k)=Int(Double(nx_p)/2.0+eluvw(0,k)*visFreq(chan)*scale(0)+0.5);
1043 0 : Int newrow=locv(chan, k)*nx_p+locu(chan, k);
1044 0 : if(rowToIndex[newrow] <0){
1045 0 : rowToIndex[newrow]=numUniq;
1046 0 : ++numUniq;
1047 : }
1048 : }
1049 : }
1050 : }
1051 : //cerr << "numrows " << vb.nRows() << " nuniq " << numUniq << endl;
1052 0 : Cube<Complex> grid(npol_p, nchan_p, numUniq);
1053 0 : Matrix<Float> wght(npol_p,numUniq);
1054 0 : Cube<Float> wghtSpec(npol_p,nchan_p,numUniq);
1055 0 : Cube<Bool> flag(npol_p, nchan_p,numUniq);
1056 0 : Vector<Bool> rowFlag(numUniq);
1057 0 : Matrix<Double> uvw(3, numUniq);
1058 0 : Vector<Int> ant1(numUniq);
1059 0 : Vector<Int> ant2(numUniq);
1060 0 : Vector<Double> timeCen(numUniq);
1061 0 : Vector<uInt> rowids(numUniq);
1062 0 : Vector<Bool> rowvisited(numUniq, false);
1063 :
1064 0 : Vector<Double> invLambda=visFreq/C::c;
1065 0 : Vector<Double> phasmult(vb.nChannels(),0.0);
1066 0 : for (Int k=0; k < nrows; ++k){
1067 0 : if(needRot)
1068 0 : phasmult=phasor(k)*invLambda;
1069 0 : for(Int chan=0; chan < vb.nChannels(); ++chan ){
1070 0 : if(chanMap_p(chan) >=0){
1071 : //Have to reject uvws out of range.
1072 0 : if(locu(chan,k) > -1 && locu(chan, k) < nx_p && locv(chan,k)> -1 && locv(chan,k) < ny_p){
1073 0 : Int newrow=locv(chan,k)*nx_p+locu(chan,k);
1074 0 : Int actrow=rowToIndex[newrow];
1075 0 : rowids[actrow]=uInt(newrow);
1076 0 : if(!rowvisited[actrow]){
1077 0 : rowvisited[actrow]=true;
1078 0 : rowFlag[actrow]=msc.flagRow()(newrow);
1079 0 : wghtSpec[actrow]=msc.weightSpectrum().get(newrow);
1080 0 : flag[actrow]=msc.flag().get(newrow);
1081 0 : uvw[actrow]=msc.uvw().get(newrow);
1082 0 : ant1[actrow]=msc.antenna1()(newrow);
1083 0 : ant2[actrow]=msc.antenna2()(newrow);
1084 0 : timeCen[actrow]=msc.time()(newrow);
1085 0 : grid[actrow]=msc.data().get(newrow);
1086 : }
1087 0 : if(rowFlag[actrow] && !(vb.flagRow()(k))){
1088 : // cerr << newrow << " rowFlag " << rowFlag[actrow] <<" rowids " << vb.rowIds()[k] << " uvw2 " << uvw(2 ,actrow) << endl;
1089 0 : rowFlag[actrow]=false;
1090 0 : uvw(2,actrow)=eluvw(2,k);
1091 : // cerr << newrow << " rowFlag " << rowFlag[actrow] <<" rowids " << vb.rowIds()[k] << " uvw2 " << uvw(2 ,actrow) << endl;
1092 0 : ant1[actrow]=vb.antenna1()(k);
1093 0 : ant2[actrow]=vb.antenna2()(k);
1094 0 : timeCen[actrow]=vb.time()(k);
1095 : }
1096 0 : for(Int pol=0; pol < vb.nCorrelations(); ++pol){
1097 0 : if((!vb.flagCube()(pol,chan, k)) && (polMap_p(pol)>=0) && (vb.weight()(pol,k)>0.0)){
1098 0 : Complex toB=hasCorrected ? vb.visCubeCorrected()(pol,chan,k)*vb.weight()(pol,k):
1099 0 : vb.visCube()(pol,chan,k)*vb.weight()(pol,k);
1100 0 : if(needRot){
1101 : Double s, c;
1102 0 : SINCOS(phasmult(chan), s, c);
1103 0 : Complex elphas(c, s);
1104 0 : toB *= elphas;
1105 : }
1106 0 : grid(polMap_p(pol),chanMap_p(chan),actrow)
1107 0 : = (grid(polMap_p(pol),chanMap_p(chan),actrow)*wghtSpec(polMap_p(pol),chanMap_p(chan),actrow)
1108 0 : + toB)/(vb.weight()(pol,k)+wghtSpec(polMap_p(pol),chanMap_p(chan),actrow));
1109 0 : flag(polMap_p(pol),chanMap_p(chan),actrow)=false;
1110 : //cerr << "weights " << max(vb.weight()) << " spec " << max(vb.weightSpectrum()) << endl;
1111 : //wghtSpec(polMap_p(pol),chanMap_p(chan), newrow)+=vb.weightSpectrum()(pol, chan, k);
1112 0 : wghtSpec(polMap_p(pol),chanMap_p(chan),actrow) += vb.weight()(pol,k);
1113 : }
1114 : }
1115 : }
1116 : }
1117 : //This should go for later
1118 : /*for (Int pol=0; pol < wght.shape()(0); ++pol){
1119 : //cerr << "shape min max "<< median(wghtSpec.xyPlane(newrow).row(pol)) << " " << min(wghtSpec.xyPlane(newrow).row(pol)) << " "<< max(wghtSpec.xyPlane(newrow).row(pol)) << endl;
1120 : wght(pol,actrow)=median(wghtSpec.xyPlane(actrow).row(pol));
1121 : //cerr << "pol " << pol << " newrow "<< newrow << " weight "<< wght(pol, newrow) << endl;
1122 : }*/
1123 : }
1124 : }
1125 : //now lets put back the stuff
1126 : //cerr << "rowids " << rowids << endl;
1127 0 : RefRows elrow(rowids);
1128 0 : msc.flagRow().putColumnCells(elrow, rowFlag);
1129 : //reference row
1130 : //Matrix<Complex>matdata(grid.xyPlane(k));
1131 0 : msc.data().putColumnCells(elrow, grid);
1132 : //Matrix<Float> matwgt(wghtSpec.xyPlane(k));
1133 0 : msc.weightSpectrum().putColumnCells(elrow, wghtSpec);
1134 : //Matrix<Bool> matflag(flag.xyPlane(k));
1135 0 : msc.flag().putColumnCells(elrow, flag);
1136 : //Vector<Double> vecuvw(uvw.column(k));
1137 0 : msc.uvw().putColumnCells(elrow, uvw);
1138 : //cerr << "ant1 " << ant1 << endl;
1139 0 : msc.antenna1().putColumnCells(elrow, ant1);
1140 0 : msc.antenna2().putColumnCells(elrow, ant2);
1141 0 : msc.time().putColumnCells(elrow, timeCen);
1142 : //outMsPtr_p->flush(true);
1143 : ///////
1144 :
1145 :
1146 0 : }
1147 :
1148 0 : void MSUVBin::inplaceSmallBW(const vi::VisBuffer2& vb){
1149 : //Dang i thought the new vb will return Data or FloatData if correctedData was
1150 : //not there
1151 0 : Bool hasCorrected=!(MSMainColumns(vb.ms()).correctedData().isNull());
1152 0 : Int nrows=vb.nRows();
1153 0 : MSColumns msc(*outMsPtr_p);
1154 0 : SpectralCoordinate spec=csys_p.spectralCoordinate(2);
1155 : Double refFreq;
1156 0 : spec.toWorld(refFreq, Double(nchan_p)/2.0);
1157 0 : DirectionCoordinate thedir=csys_p.directionCoordinate(0);
1158 0 : Vector<Float> scale(2);
1159 0 : scale(0)=(nx_p*thedir.increment()(0))/C::c;
1160 0 : scale(1)=(ny_p*thedir.increment()(1))/C::c;
1161 : ///Index
1162 0 : Vector<Int> rowToIndex(nx_p*ny_p, -1);
1163 0 : Vector<Int> locu(nrows, -1);
1164 0 : Vector<Int> locv(nrows, -1);
1165 0 : Vector<Double> visFreq=vb.getFrequencies(0, MFrequency::LSRK);
1166 0 : Int numUniq=0;
1167 0 : Vector<Double> phasor;
1168 0 : Matrix<Double> eluvw;
1169 0 : eluvw=vb.uvw();
1170 0 : Bool needRot=vbutil_p.rotateUVW(vb, phaseCenter_p, eluvw, phasor);
1171 : //cerr << "SHAPES " << eluvw.shape() << " " << vb.uvw().shape() << endl;
1172 0 : for (Int k=0; k <nrows; ++k){
1173 0 : locv(k)=Int(Double(ny_p)/2.0+eluvw(1,k)*refFreq*scale(1)+0.5);
1174 0 : locu(k)=Int(Double(nx_p)/2.0+eluvw(0,k)*refFreq*scale(0)+0.5);
1175 : //cerr << "locu " << locu(k) << " locv " << locv(k) << " eluvw " << eluvw(0,k) << " uvw " << vb.uvw()(0,k) << endl;
1176 0 : if(locu(k) > -1 && locu(k) < nx_p && locv(k)> -1 && locv(k) < ny_p){
1177 0 : Int newrow=locv(k)*nx_p+locu(k);
1178 0 : if(rowToIndex[newrow] <0){
1179 0 : rowToIndex[newrow]=numUniq;
1180 0 : ++numUniq;
1181 : }
1182 : }
1183 :
1184 : }
1185 :
1186 : //cerr << "numrows " << vb.nRows() << " nuniq " << numUniq << endl;
1187 0 : Cube<Complex> grid(npol_p, nchan_p, numUniq);
1188 0 : Matrix<Float> wght(npol_p,numUniq);
1189 0 : Cube<Float> wghtSpec(npol_p,nchan_p,numUniq);
1190 0 : Cube<Bool> flag(npol_p, nchan_p,numUniq);
1191 0 : Vector<Bool> rowFlag(numUniq);
1192 0 : Matrix<Double> uvw(3, numUniq);
1193 0 : Vector<Int> ant1(numUniq);
1194 0 : Vector<Int> ant2(numUniq);
1195 0 : Vector<Double> timeCen(numUniq);
1196 0 : Vector<uInt> rowids(numUniq);
1197 0 : Vector<Bool> rowvisited(numUniq, false);
1198 :
1199 0 : Vector<Double> invLambda=visFreq/C::c;
1200 0 : for (Int k=0; k < nrows; ++k){
1201 :
1202 0 : if(locu(k) > -1 && locu(k) < nx_p && locv(k)> -1 && locv(k) < ny_p){
1203 0 : Int newrow=locv(k)*nx_p+locu(k);
1204 0 : Int actrow=rowToIndex[newrow];
1205 0 : rowids[actrow]=uInt(newrow);
1206 0 : if(!rowvisited[actrow]){
1207 0 : rowvisited[actrow]=true;
1208 0 : rowFlag[actrow]=msc.flagRow()(newrow);
1209 0 : wghtSpec[actrow]=msc.weightSpectrum().get(newrow);
1210 0 : flag[actrow]=msc.flag().get(newrow);
1211 0 : uvw[actrow]=msc.uvw().get(newrow);
1212 0 : ant1[actrow]=msc.antenna1()(newrow);
1213 0 : ant2[actrow]=msc.antenna2()(newrow);
1214 0 : timeCen[actrow]=msc.time()(newrow);
1215 0 : grid[actrow]=msc.data().get(newrow);
1216 : }
1217 0 : if(rowFlag[actrow] && !(vb.flagRow()(k))){
1218 : // cerr << newrow << " rowFlag " << rowFlag[actrow] <<" rowids " << vb.rowIds()[k] << " uvw2 " << uvw(2 ,actrow) << endl;
1219 0 : rowFlag[actrow]=false;
1220 0 : uvw(2,actrow)=eluvw(2,k);
1221 : // cerr << newrow << " rowFlag " << rowFlag[actrow] <<" rowids " << vb.rowIds()[k] << " uvw2 " << uvw(2 ,actrow) << endl;
1222 0 : ant1[actrow]=vb.antenna1()(k);
1223 0 : ant2[actrow]=vb.antenna2()(k);
1224 0 : timeCen[actrow]=vb.time()(k);
1225 : }
1226 0 : Complex elphas(1.0, 0.0);
1227 0 : for(Int chan=0; chan < vb.nChannels(); ++chan ){
1228 :
1229 0 : if(chanMap_p(chan) >=0){
1230 0 : if(needRot){
1231 0 : Double phasmult=phasor(k)*invLambda(chan);
1232 : Double s, c;
1233 0 : SINCOS(phasmult, s, c);
1234 0 : elphas=Complex(c, s);
1235 : }
1236 :
1237 0 : for(Int pol=0; pol < vb.nCorrelations(); ++pol){
1238 0 : if((!vb.flagCube()(pol,chan, k)) && (polMap_p(pol)>=0) && (vb.weight()(pol,k)>0.0)){
1239 0 : Complex toB=hasCorrected ? vb.visCubeCorrected()(pol,chan,k)*vb.weight()(pol,k):
1240 0 : vb.visCube()(pol,chan,k)*vb.weight()(pol,k);
1241 0 : if(needRot)
1242 0 : toB *= elphas;
1243 0 : grid(polMap_p(pol),chanMap_p(chan),actrow)
1244 0 : = (grid(polMap_p(pol),chanMap_p(chan),actrow)*wghtSpec(polMap_p(pol),chanMap_p(chan),actrow)
1245 0 : + toB)/(vb.weight()(pol,k)+wghtSpec(polMap_p(pol),chanMap_p(chan),actrow));
1246 0 : flag(polMap_p(pol),chanMap_p(chan),actrow)=false;
1247 : //cerr << "weights " << max(vb.weight()) << " spec " << max(vb.weightSpectrum()) << endl;
1248 : //wghtSpec(polMap_p(pol),chanMap_p(chan), newrow)+=vb.weightSpectrum()(pol, chan, k);
1249 0 : wghtSpec(polMap_p(pol),chanMap_p(chan),actrow) += vb.weight()(pol,k);
1250 : }
1251 : }
1252 : }
1253 : }
1254 : //This should go for later
1255 : /*for (Int pol=0; pol < wght.shape()(0); ++pol){
1256 : //cerr << "shape min max "<< median(wghtSpec.xyPlane(newrow).row(pol)) << " " << min(wghtSpec.xyPlane(newrow).row(pol)) << " "<< max(wghtSpec.xyPlane(newrow).row(pol)) << endl;
1257 : wght(pol,actrow)=median(wghtSpec.xyPlane(actrow).row(pol));
1258 : //cerr << "pol " << pol << " newrow "<< newrow << " weight "<< wght(pol, newrow) << endl;
1259 : }*/
1260 : }
1261 : }
1262 : //now lets put back the stuff
1263 0 : RefRows elrow(rowids);
1264 : //cerr << "ROWIDS " << rowids << endl;
1265 0 : msc.flagRow().putColumnCells(elrow, rowFlag);
1266 : //reference row
1267 : //Matrix<Complex>matdata(grid.xyPlane(k));
1268 0 : msc.data().putColumnCells(elrow, grid);
1269 : //Matrix<Float> matwgt(wghtSpec.xyPlane(k));
1270 0 : msc.weightSpectrum().putColumnCells(elrow, wghtSpec);
1271 : //Matrix<Bool> matflag(flag.xyPlane(k));
1272 0 : msc.flag().putColumnCells(elrow, flag);
1273 : //Vector<Double> vecuvw(uvw.column(k));
1274 0 : msc.uvw().putColumnCells(elrow, uvw);
1275 : //cerr << "ant1 " << ant1 << endl;
1276 0 : msc.antenna1().putColumnCells(elrow, ant1);
1277 0 : msc.antenna2().putColumnCells(elrow, ant2);
1278 0 : msc.time().putColumnCells(elrow, timeCen);
1279 : //outMsPtr_p->flush(true);
1280 : ///////
1281 :
1282 :
1283 0 : }
1284 :
1285 :
1286 :
1287 0 : void MSUVBin::gridData(const vi::VisBuffer2& vb, Cube<Complex>& grid,
1288 : Matrix<Float>& /*wght*/, Cube<Float>& wghtSpec,
1289 : Cube<Bool>& flag, Vector<Bool>& rowFlag, Matrix<Double>& uvw, Vector<Int>& ant1,
1290 : Vector<Int>& ant2, Vector<Double>& timeCen, const Matrix<Int>& /*locuv*/){
1291 : //all pixel that are touched the flag and flag Row shall be unset and the w be assigned
1292 : //later we'll deal with multiple w for the same uv
1293 : //we need polmap and chanmap;
1294 :
1295 : Double fracbw;
1296 0 : if(!datadescMap(vb, fracbw)) return;
1297 : //cerr << "fracbw " << fracbw << endl;
1298 0 : SpectralCoordinate spec=csys_p.spectralCoordinate(2);
1299 0 : DirectionCoordinate thedir=csys_p.directionCoordinate(0);
1300 0 : Double refFreq=SpectralImageUtil::worldFreq(csys_p, Double(nchan_p/2));
1301 : //Double refFreq=SpectralImageUtil::worldFreq(csys_p, Double(0));
1302 0 : Vector<Float> scale(2);
1303 0 : scale(0)=fabs(nx_p*thedir.increment()(0))/C::c;
1304 0 : scale(1)=fabs(ny_p*thedir.increment()(1))/C::c;
1305 : //Dang i thought the new vb will return Data or FloatData if correctedData was
1306 : //not there
1307 0 : Bool hasCorrected=!(MSMainColumns(vb.ms()).correctedData().isNull());
1308 : //locateuvw(locuv, vb.uvw());
1309 0 : Vector<Double> visFreq=vb.getFrequencies(0, MFrequency::LSRK);
1310 0 : for (rownr_t k=0; k < vb.nRows(); ++k){
1311 0 : if(!vb.flagRow()[k]){
1312 : Int locu, locv;
1313 : {
1314 0 : locv=Int(Double(ny_p)/2.0+vb.uvw()(1,k)*refFreq*scale(1)+0.5);
1315 0 : locu=Int(Double(nx_p)/2.0+vb.uvw()(0,k)*refFreq*scale(0)+0.5);
1316 :
1317 : }
1318 : //cerr << "fracbw " << fracbw << " pixel u " << Double(locu-nx_p/2)/refFreq/scale(0) << " data u" << vb.uvw()(0,k) <<
1319 : // " pixel v "<< Double(locv-ny_p/2)/refFreq/scale(1) << " data v "<< vb.uvw()(1,k) << endl;
1320 : //Double newU= Double(locu-nx_p/2)/refFreq/scale(0);
1321 : //Double newV= Double(locv-ny_p/2)/refFreq/scale(1);
1322 : //Double phaseCorr=((newU/vb.uvw()(0,k)-1)+(newV/vb.uvw()(1,k)-1))*vb.uvw()(2,k)*2.0*C::pi*refFreq/C::c;
1323 : //Double phaseCorr=((newU-vb.uvw()(0,k))+(newV-vb.uvw()(1,k)))*2.0*C::pi*refFreq/C::c;
1324 : //cerr << "fracbw " << fracbw << " pixel u " << Double(locu-nx_p/2)/refFreq/scale(0) << " data u" << vb.uvw()(0,k) <<
1325 : // " pixel v "<< Double(locv-ny_p/2)/refFreq/scale(1) << " data v "<< vb.uvw()(1,k) << " phase " << phaseCorr << endl;
1326 0 : for(Int chan=0; chan < vb.nChannels(); ++chan ){
1327 0 : if(chanMap_p(chan) >=0){
1328 : //Double outChanFreq;
1329 : //spec.toWorld(outChanFreq, Double(chanMap_p(chan)));
1330 0 : if(fracbw > 0.05)
1331 : {
1332 0 : locv=Int(Double(ny_p)/2.0+vb.uvw()(1,k)*visFreq(chan)*scale(1)+0.5);
1333 0 : locu=Int(Double(nx_p)/2.0+vb.uvw()(0,k)*visFreq(chan)*scale(0)+0.5);
1334 : }
1335 0 : if(locv < ny_p && locu < nx_p){
1336 0 : Int newrow=locv*nx_p+locu;
1337 0 : if(rowFlag(newrow) && !(vb.flagRow()(k))){
1338 0 : rowFlag(newrow)=false;
1339 : /////TEST
1340 : //uvw(0,newrow)=vb.uvw()(0,k);
1341 : //uvw(1,newrow)=vb.uvw()(1,k);
1342 : /////
1343 0 : uvw(2,newrow)=vb.uvw()(2,k);
1344 : //cerr << newrow << " rowids " << vb.rowIds()[k] << " uvw2 " << uvw(2 ,newrow) << endl;
1345 0 : ant1(newrow)=vb.antenna1()(k);
1346 0 : ant2(newrow)=vb.antenna2()(k);
1347 0 : timeCen(newrow)=vb.time()(k);
1348 :
1349 : }
1350 0 : for(Int pol=0; pol < vb.nCorrelations(); ++pol){
1351 0 : if((!vb.flagCube()(pol,chan, k)) && (polMap_p(pol)>=0) && (vb.weight()(pol,k)>0.0)){
1352 : // Double newU= Double(locu-nx_p/2)/refFreq/scale(0);
1353 : // Double newV= Double(locv-ny_p/2)/refFreq/scale(1);
1354 : // Double phaseCorr=((newU/vb.uvw()(0,k)-1)+(newV/vb.uvw()(1,k)-1))*vb.uvw()(2,k)*2.0*C::pi*refFreq/C::c;
1355 0 : Complex toB=hasCorrected ? vb.visCubeCorrected()(pol,chan,k)*vb.weight()(pol,k):
1356 0 : vb.visCube()(pol,chan,k)*vb.weight()(pol,k);
1357 : // Double s, c;
1358 : //SINCOS(phaseCorr, s, c);
1359 : // toB=toB*Complex(c,s);
1360 0 : grid(polMap_p(pol),chanMap_p(chan), newrow)
1361 0 : = (grid(polMap_p(pol),chanMap_p(chan), newrow)*wghtSpec(polMap_p(pol),chanMap_p(chan),newrow)
1362 0 : + toB)/(vb.weight()(pol,k)+wghtSpec(polMap_p(pol),chanMap_p(chan),newrow));
1363 0 : flag(polMap_p(pol),chanMap_p(chan), newrow)=false;
1364 : //cerr << "weights " << max(vb.weight()) << " spec " << max(vb.weightSpectrum()) << endl;
1365 : //wghtSpec(polMap_p(pol),chanMap_p(chan), newrow)+=vb.weightSpectrum()(pol, chan, k);
1366 0 : wghtSpec(polMap_p(pol),chanMap_p(chan), newrow) += vb.weight()(pol,k);
1367 : }
1368 : ///We should do that at the end totally
1369 : //wght(pol,newrow)=median(wghtSpec.xyPlane(newrow).row(pol));
1370 : }
1371 : }//locu && locv
1372 : }
1373 : }
1374 : //sum wgtspec along channels for weight
1375 : /*for (Int pol=0; pol < wght.shape()(0); ++pol){
1376 : //cerr << "shape min max "<< median(wghtSpec.xyPlane(newrow).row(pol)) << " " << min(wghtSpec.xyPlane(newrow).row(pol)) << " "<< max(wghtSpec.xyPlane(newrow).row(pol)) << endl;
1377 : wght(pol,newrow)=median(wghtSpec.xyPlane(newrow).row(pol));
1378 : //cerr << "pol " << pol << " newrow "<< newrow << " weight "<< wght(pol, newrow) << endl;
1379 : }
1380 : */
1381 : }
1382 : }
1383 :
1384 :
1385 :
1386 :
1387 0 : }
1388 :
1389 0 : void MSUVBin::gridData(const vi::VisBuffer2& vb, Cube<Complex>& grid,
1390 : Matrix<Float>& /*wght*/, Cube<Float>& wghtSpec,
1391 : Cube<Bool>& flag, Vector<Bool>& rowFlag, Matrix<Double>& uvw, Vector<Int>& ant1,
1392 : Vector<Int>& ant2, Vector<Double>& timeCen, const Int startchan, const Int endchan){
1393 : //all pixel that are touched the flag and flag Row shall be unset and the w be assigned
1394 : //later we'll deal with multiple w for the same uv
1395 : //we need polmap and chanmap;
1396 :
1397 : Double fracbw;
1398 0 : if(!datadescMap(vb, fracbw)) return;
1399 : //cerr << "fracbw " << fracbw << endl;
1400 0 : SpectralCoordinate spec=csys_p.spectralCoordinate(2);
1401 0 : DirectionCoordinate thedir=csys_p.directionCoordinate(0);
1402 0 : Double refFreq=SpectralImageUtil::worldFreq(csys_p, Double(nchan_p/2));
1403 : //Double refFreq=SpectralImageUtil::worldFreq(csys_p, Double(0));
1404 0 : Vector<Float> scale(2);
1405 0 : scale(0)=fabs(Double(nx_p)*thedir.increment()(0))/C::c;
1406 0 : scale(1)=fabs(Double(ny_p)*thedir.increment()(1))/C::c;
1407 : //cerr << "dir incr "<< thedir.increment() << " nx ny " << nx_p << " " << ny_p << endl;
1408 : //cerr << "SCALE " << scale << " fracbw " << fracbw << endl;
1409 :
1410 : //cerr << "chanmap " << chanMap_p << " pol map " <<polMap_p << " weight " << wghtSpec.shape() << endl;
1411 : //Dang i thought the new vb will return Data or FloatData if correctedData was
1412 : //not there
1413 0 : Bool hasCorrected=!(MSMainColumns(vb.ms()).correctedData().isNull());
1414 : //locateuvw(locuv, vb.uvw());
1415 0 : Vector<Double> visFreq=vb.getFrequencies(0, MFrequency::LSRK);
1416 0 : for (rownr_t k=0; k < vb.nRows(); ++k){
1417 0 : if(!vb.flagRow()[k]){
1418 : Int locu, locv;
1419 : {
1420 0 : locv=Int(Double(ny_p)/2.0+vb.uvw()(1,k)*refFreq*scale(1)+0.5);
1421 0 : locu=Int(Double(nx_p)/2.0+vb.uvw()(0,k)*refFreq*scale(0)+0.5);
1422 :
1423 : }
1424 : //cerr << "fracbw " << fracbw << " pixel u " << Double(locu-nx_p/2)/refFreq/scale(0) << " data u" << vb.uvw()(0,k) <<
1425 : // " pixel v "<< Double(locv-ny_p/2)/refFreq/scale(1) << " data v "<< vb.uvw()(1,k) << endl;
1426 : //Double newU= Double(locu-nx_p/2)/refFreq/scale(0);
1427 : //Double newV= Double(locv-ny_p/2)/refFreq/scale(1);
1428 : //Double phaseCorr=((newU/vb.uvw()(0,k)-1)+(newV/vb.uvw()(1,k)-1))*vb.uvw()(2,k)*2.0*C::pi*refFreq/C::c;
1429 : //Double phaseCorr=((newU-vb.uvw()(0,k))+(newV-vb.uvw()(1,k)))*2.0*C::pi*refFreq/C::c;
1430 : //cerr << "fracbw " << fracbw << " pixel u " << Double(locu-nx_p/2)/refFreq/scale(0) << " data u" << vb.uvw()(0,k) <<
1431 : // " pixel v "<< Double(locv-ny_p/2)/refFreq/scale(1) << " data v "<< vb.uvw()(1,k) << " phase " << phaseCorr << endl;
1432 0 : for(Int chan=0; chan < vb.nChannels(); ++chan ){
1433 0 : if(chanMap_p(chan) >=startchan && chanMap_p(chan) <=endchan){
1434 : //Double outChanFreq;
1435 : //spec.toWorld(outChanFreq, Double(chanMap_p(chan)));
1436 0 : if(fracbw > 0.05)
1437 : {
1438 0 : locv=Int(Double(ny_p)/2.0+vb.uvw()(1,k)*visFreq(chan)*scale(1)+0.5);
1439 0 : locu=Int(Double(nx_p)/2.0+vb.uvw()(0,k)*visFreq(chan)*scale(0)+0.5);
1440 : }
1441 0 : if(locv < ny_p && locu < nx_p && locv >=0 && locu >=0){
1442 0 : Int newrow=locv*nx_p+locu;
1443 0 : if(rowFlag(newrow) && !(vb.flagRow()(k))){
1444 0 : rowFlag(newrow)=false;
1445 : /////TEST
1446 : //uvw(0,newrow)=vb.uvw()(0,k);
1447 : //uvw(1,newrow)=vb.uvw()(1,k);
1448 : /////
1449 0 : uvw(2,newrow)=vb.uvw()(2,k);
1450 : //cerr << newrow << " rowids " << vb.rowIds()[k] << " uvw2 " << uvw(2 ,newrow) << endl;
1451 0 : ant1(newrow)=vb.antenna1()(k);
1452 0 : ant2(newrow)=vb.antenna2()(k);
1453 0 : timeCen(newrow)=vb.time()(k);
1454 :
1455 : }
1456 0 : for(Int pol=0; pol < vb.nCorrelations(); ++pol){
1457 0 : if((!vb.flagCube()(pol,chan, k)) && (polMap_p(pol)>=0) && (vb.weight()(pol,k)>0.0)){
1458 : // Double newU= Double(locu-nx_p/2)/refFreq/scale(0);
1459 : // Double newV= Double(locv-ny_p/2)/refFreq/scale(1);
1460 : // Double phaseCorr=((newU/vb.uvw()(0,k)-1)+(newV/vb.uvw()(1,k)-1))*vb.uvw()(2,k)*2.0*C::pi*refFreq/C::c;
1461 0 : Complex toB=hasCorrected ? vb.visCubeCorrected()(pol,chan,k)*vb.weight()(pol,k):
1462 0 : vb.visCube()(pol,chan,k)*vb.weight()(pol,k);
1463 : // Double s, c;
1464 : //SINCOS(phaseCorr, s, c);
1465 : // toB=toB*Complex(c,s);
1466 0 : grid(polMap_p(pol),chanMap_p(chan)-startchan, newrow)
1467 0 : = (grid(polMap_p(pol),chanMap_p(chan)-startchan, newrow)
1468 0 : + toB); ///(vb.weight()(pol,k)+wghtSpec(polMap_p(pol),chanMap_p(chan)-startchan,newrow));
1469 0 : flag(polMap_p(pol),chanMap_p(chan)-startchan, newrow)=false;
1470 : //cerr << "weights " << max(vb.weight()) << " spec " << max(vb.weightSpectrum()) << endl;
1471 : //wghtSpec(polMap_p(pol),chanMap_p(chan), newrow)+=vb.weightSpectrum()(pol, chan, k);
1472 0 : wghtSpec(polMap_p(pol),chanMap_p(chan)-startchan, newrow) += vb.weight()(pol,k);
1473 : }
1474 : ///We should do that at the end totally
1475 : //wght(pol,newrow)=median(wghtSpec.xyPlane(newrow).row(pol));
1476 : }
1477 : }//locu && locv
1478 : }
1479 : }
1480 : //sum wgtspec along channels for weight
1481 : /*for (Int pol=0; pol < wght.shape()(0); ++pol){
1482 : //cerr << "shape min max "<< median(wghtSpec.xyPlane(newrow).row(pol)) << " " << min(wghtSpec.xyPlane(newrow).row(pol)) << " "<< max(wghtSpec.xyPlane(newrow).row(pol)) << endl;
1483 : wght(pol,newrow)=median(wghtSpec.xyPlane(newrow).row(pol));
1484 : //cerr << "pol " << pol << " newrow "<< newrow << " weight "<< wght(pol, newrow) << endl;
1485 : }
1486 : */
1487 : }
1488 : }
1489 :
1490 :
1491 :
1492 :
1493 0 : }
1494 0 : void MSUVBin::gridDataConv(const vi::VisBuffer2& vb, Cube<Complex>& grid,
1495 : Matrix<Float>& /*wght*/, Cube<Complex>& wghtSpec,
1496 : Cube<Bool>& flag, Vector<Bool>& rowFlag, Matrix<Double>& uvw, Vector<Int>& ant1,
1497 : Vector<Int>& ant2, Vector<Double>& timeCen, const Int startchan, const Int endchan, const Cube<Complex>& convFunc, const Vector<Int>& convSupport, const Double wScale, const Int convSampling){
1498 : //all pixel that are touched the flag and flag Row shall be unset and the w be assigned
1499 : //later we'll deal with multiple w for the same uv
1500 : //we need polmap and chanmap;
1501 :
1502 : Double fracbw;
1503 0 : if(!datadescMap(vb, fracbw)) return;
1504 : //cerr << "fracbw " << fracbw << endl;
1505 0 : SpectralCoordinate spec=csys_p.spectralCoordinate(2);
1506 0 : DirectionCoordinate thedir=csys_p.directionCoordinate(0);
1507 0 : Double refFreq=SpectralImageUtil::worldFreq(csys_p, Double(nchan_p/2));
1508 : //Double refFreq=SpectralImageUtil::worldFreq(csys_p, Double(0));
1509 0 : Vector<Float> scale(2);
1510 0 : scale(0)=fabs(Double(nx_p)*thedir.increment()(0))/C::c;
1511 0 : scale(1)=fabs(Double(ny_p)*thedir.increment()(1))/C::c;
1512 : //Dang i thought the new vb will return Data or FloatData if correctedData was
1513 : //not there
1514 0 : Bool hasCorrected=!(MSMainColumns(vb.ms()).correctedData().isNull());
1515 : //locateuvw(locu v, vb.uvw());
1516 0 : Vector<Double> visFreq=vb.getFrequencies(0, MFrequency::LSRK);
1517 : //cerr << "support " << convSupport << endl;
1518 0 : Vector<Double> phasor;
1519 0 : Matrix<Double> eluvw;
1520 0 : eluvw=vb.uvw();
1521 0 : Bool needRot=vbutil_p.rotateUVW(vb, phaseCenter_p, eluvw, phasor);
1522 0 : Vector<Double> invLambda=visFreq/C::c;
1523 0 : for (rownr_t k=0; k < vb.nRows(); ++k){
1524 0 : if(!vb.flagRow()[k]){
1525 : Int locu, locv, locw,supp, offu, offv;
1526 : {
1527 0 : locv=Int(Double(ny_p)/2.0+vb.uvw()(1,k)*refFreq*scale(1)+0.5);
1528 0 : locu=Int(Double(nx_p)/2.0+vb.uvw()(0,k)*refFreq*scale(0)+0.5);
1529 0 : offv=Int ((Double(locv)-(Double(ny_p)/2.0+vb.uvw()(1,k)*refFreq*scale(1)))*Double(convSampling)+0.5);
1530 0 : offu=Int ((Double(locu)-(Double(nx_p)/2.0+vb.uvw()(0,k)*refFreq*scale(0)))*Double(convSampling)+0.5);
1531 0 : locw=Int(sqrt(fabs(wScale*vb.uvw()(2,k)*refFreq/C::c))+0.5);
1532 : //cerr << "locw " << locw << " " << " offs " << offu << " " << offv << endl;
1533 0 : supp=locw < convSupport.shape()[0] ? convSupport(locw) :convSupport(convSupport.nelements()-1) ;
1534 : //////////////////
1535 : //supp=10;
1536 : /////////////////
1537 : }
1538 : //cerr << "fracbw " << fracbw << " pixel u " << Double(locu-nx_p/2)/refFreq/scale(0) << " data u" << vb.uvw()(0,k) <<
1539 : // " pixel v "<< Double(locv-ny_p/2)/refFreq/scale(1) << " data v "<< vb.uvw()(1,k) << endl;
1540 : //Double newU= Double(locu-nx_p/2)/refFreq/scale(0);
1541 : //Double newV= Double(locv-ny_p/2)/refFreq/scale(1);
1542 : //Double phaseCorr=((newU/vb.uvw()(0,k)-1)+(newV/vb.uvw()(1,k)-1))*vb.uvw()(2,k)*2.0*C::pi*refFreq/C::c;
1543 : //Double phaseCorr=((newU-vb.uvw()(0,k))+(newV-vb.uvw()(1,k)))*2.0*C::pi*refFreq/C::c;
1544 : //cerr << "fracbw " << fracbw << " pixel u " << Double(locu-nx_p/2)/refFreq/scale(0) << " data u" << vb.uvw()(0,k) <<
1545 : // " pixel v "<< Double(locv-ny_p/2)/refFreq/scale(1) << " data v "<< vb.uvw()(1,k) << " phase " << phaseCorr << endl;
1546 0 : Vector<Int> newrow((2*supp+1)*(2*supp+1), -1);
1547 0 : if(locv < ny_p && locu < nx_p){
1548 0 : for (Int yy=0; yy< Int(2*supp+1); ++yy){
1549 0 : Int newlocv=locv+ yy-supp;
1550 0 : if(newlocv >0 && newlocv < ny_p){
1551 0 : for (Int xx=0; xx< Int(2*supp+1); ++xx){
1552 0 : Int newlocu=locu+xx-supp;
1553 0 : if(newlocu < nx_p && newlocu >0){
1554 0 : newrow(yy*(2*supp+1)+xx)=newlocv*nx_p+newlocu;
1555 : }
1556 : }
1557 : }
1558 : }
1559 : }
1560 0 : Complex elphas(1.0, 0.0);
1561 0 : for(Int chan=0; chan < vb.nChannels(); ++chan ){
1562 0 : if(chanMap_p(chan) >=startchan && chanMap_p(chan) <=endchan){
1563 : //Double outChanFreq;
1564 : //spec.toWorld(outChanFreq, Double(chanMap_p(chan)));
1565 0 : if(needRot){
1566 0 : Double phasmult=phasor(k)*invLambda(chan);
1567 : Double s, c;
1568 0 : SINCOS(phasmult, s, c);
1569 0 : elphas=Complex(c, s);
1570 : }
1571 0 : if(fracbw > 0.05)
1572 : {
1573 0 : locv=Int(Double(ny_p)/2.0+vb.uvw()(1,k)*visFreq(chan)*scale(1)+0.5);
1574 0 : locu=Int(Double(nx_p)/2.0+vb.uvw()(0,k)*visFreq(chan)*scale(0)+0.5);
1575 : }
1576 0 : if(locv < ny_p && locu < nx_p){
1577 0 : for (Int yy=0; yy< Int(2*supp+1); ++yy){
1578 0 : Int locy=abs((yy-supp)*convSampling+offv);
1579 : //cerr << "y " << yy << " locy " << locy ;
1580 0 : for (Int xx=0; xx< Int(2*supp+1); ++xx){
1581 0 : Int locx=abs((xx-supp)*convSampling+offu);
1582 : //cerr << " x " << xx << " locx " << locx << endl;
1583 0 : Int jj=yy*(2*supp+1)+xx;
1584 0 : if(newrow[jj] >-1){
1585 0 : Complex cwt=convFunc(locx, locy, locw);
1586 : ////////////////TEST
1587 : //cwt=Complex(1.0, 0.0);
1588 : ///////////////////////////////////
1589 0 : if(vb.uvw()(2,k) > 0.0)
1590 0 : cwt=conj(cwt);
1591 0 : if(rowFlag(newrow[jj]) && !(vb.flagRow()(k))){
1592 0 : rowFlag(newrow[jj])=false;
1593 : /////TEST
1594 : //uvw(0,newrow)=vb.uvw()(0,k);
1595 : //uvw(1,newrow)=vb.uvw()(1,k);
1596 : /////
1597 0 : uvw(2,newrow[jj])=0;
1598 : //cerr << newrow << " rowids " << vb.rowIds()[k] << " uvw2 " << uvw(2 ,newrow) << endl;
1599 0 : ant1(newrow[jj])=vb.antenna1()(k);
1600 0 : ant2(newrow[jj])=vb.antenna2()(k);
1601 0 : timeCen(newrow[jj])=vb.time()(k);
1602 :
1603 : }
1604 0 : Int lechan=chanMap_p(chan)-startchan;
1605 0 : for(Int pol=0; pol < vb.nCorrelations(); ++pol){
1606 0 : if((!vb.flagCube()(pol,chan, k)) && (polMap_p(pol)>=0) && (vb.weight()(pol,k)>0.0) && (fabs(cwt) > 0.0)){
1607 : // Double newU= Double(locu-nx_p/2)/refFreq/scale(0);
1608 : // Double newV= Double(locv-ny_p/2)/refFreq/scale(1);
1609 : // Double phaseCorr=((newU/vb.uvw()(0,k)-1)+(newV/vb.uvw()(1,k)-1))*vb.uvw()(2,k)*2.0*C::pi*refFreq/C::c;
1610 : Complex toB=hasCorrected ?
1611 0 : vb.visCubeCorrected()(pol,chan,k)*vb.weight()(pol,k)*cwt:
1612 0 : vb.visCube()(pol,chan,k)*vb.weight()(pol,k)*cwt;
1613 : //Complex dat=hasCorrected ? vb.visCubeCorrected()(pol,chan,k) :
1614 : // vb.visCube()(pol,chan,k);
1615 0 : if(needRot)
1616 0 : toB *=elphas;
1617 : // Double s, c;
1618 : //SINCOS(phaseCorr, s, c);
1619 : // toB=toB*Complex(c,s);
1620 : //Float elwgt=vb.weight()(pol,k)* fabs(real(cwt));
1621 : //////////////TESTING
1622 0 : Complex elwgt=vb.weight()(pol,k)* cwt;
1623 : //Complex elwgt=real(dat) !=0.0 ? real(toB)/real(dat): 0.0;
1624 : ///////////////////////////
1625 0 : grid(polMap_p(pol),lechan, newrow[jj])
1626 0 : = (grid(polMap_p(pol),lechan, newrow[jj])/*wghtSpec(polMap_p(pol),lechan,newrow[jj])*/
1627 0 : + toB);///(elwgt+wghtSpec(polMap_p(pol),lechan,newrow[jj]));
1628 0 : flag(polMap_p(pol), lechan, newrow[jj])=false;
1629 : //cerr << "weights " << max(vb.weight()) << " spec " << max(vb.weightSpectrum()) << endl;
1630 : //wghtSpec(polMap_p(pol),chanMap_p(chan), newrow)+=vb.weightSpectrum()(pol, chan, k);
1631 : // if( wghtSpec(polMap_p(pol),lechan, newrow[jj]) > 10)
1632 : //if(newrow[jj]==2149026)
1633 : // cerr <<xx << " : " << yy << " : " <<jj <<" spec " << wghtSpec(polMap_p(pol),lechan, newrow[jj]) << " elwgt " << elwgt << " chan " <<chan << " : " << lechan << " pol " << pol << " newrow " << newrow[jj] << " vbrow " << vb.rowIds()(k) << " uv " << uvw(0, newrow[jj]) << " " << uvw(1, newrow[jj]) << " " << vb.uvw()(0,k) << " " << vb.uvw()(1,k) <<endl;
1634 0 : wghtSpec(polMap_p(pol),lechan, newrow[jj]) += elwgt;
1635 : }
1636 : ///We should do that at the end totally
1637 : //wght(pol,newrow)=median(wghtSpec.xyPlane(newrow).row(pol));
1638 : }
1639 : }//if(newrow)
1640 : }
1641 : }
1642 : }//locu && locv
1643 : }
1644 : }
1645 : //sum wgtspec along channels for weight
1646 : /*for (Int pol=0; pol < wght.shape()(0); ++pol){
1647 : //cerr << "shape min max "<< median(wghtSpec.xyPlane(newrow).row(pol)) << " " << min(wghtSpec.xyPlane(newrow).row(pol)) << " "<< max(wghtSpec.xyPlane(newrow).row(pol)) << endl;
1648 : wght(pol,newrow)=median(wghtSpec.xyPlane(newrow).row(pol));
1649 : //cerr << "pol " << pol << " newrow "<< newrow << " weight "<< wght(pol, newrow) << endl;
1650 : }
1651 : */
1652 0 : }
1653 : }
1654 :
1655 :
1656 :
1657 :
1658 0 : }
1659 0 : void MSUVBin::gridDataConvThr(const vi::VisBuffer2& vb, Cube<Complex>& grid,
1660 : Cube<Complex>& wghtSpec,
1661 : Cube<Bool>& flag, Vector<Bool>& rowFlag, Matrix<Double>& uvw, Vector<Int>& ant1,
1662 : Vector<Int>& ant2, Vector<Double>& timeCen, const Int startchan, const Int endchan, const Cube<Complex>& convFunc, const Vector<Int>& convSupport, const Double wScale, const Int convSampling){
1663 : //all pixel that are touched the flag and flag Row shall be unset and the w be assigned
1664 : //later we'll deal with multiple w for the same uv
1665 : //we need polmap and chanmap;
1666 :
1667 : Double fracbw;
1668 0 : if(!datadescMap(vb, fracbw)) return;
1669 : //cerr << "fracbw " << fracbw << endl;
1670 0 : SpectralCoordinate spec=csys_p.spectralCoordinate(2);
1671 0 : DirectionCoordinate thedir=csys_p.directionCoordinate(0);
1672 0 : Double refFreq=SpectralImageUtil::worldFreq(csys_p, Double(nchan_p/2));
1673 : //Double refFreq=SpectralImageUtil::worldFreq(csys_p, Double(0));
1674 0 : Vector<Float> scale(2);
1675 0 : scale(0)=fabs(Double(nx_p)*thedir.increment()(0))/C::c;
1676 0 : scale(1)=fabs(Double(ny_p)*thedir.increment()(1))/C::c;
1677 : //Dang i thought the new vb will return Data or FloatData if correctedData was
1678 : //not there
1679 0 : Bool hasCorrected=!(MSMainColumns(vb.ms()).correctedData().isNull());
1680 0 : Vector<Double> visFreq=vb.getFrequencies(0, MFrequency::LSRK);
1681 : //cerr << "support " << convSupport << endl;
1682 0 : Vector<Double> phasor;
1683 0 : Matrix<Double> eluvw;
1684 0 : eluvw=vb.uvw();
1685 0 : Bool needRot=vbutil_p.rotateUVW(vb, phaseCenter_p, eluvw, phasor);
1686 : Bool gridCopy, weightCopy, flagCopy, rowFlagCopy, uvwCopy, ant1Copy, ant2Copy, timeCenCopy, sumWeightCopy, numVisCopy;
1687 0 : Complex * gridStor=grid.getStorage(gridCopy);
1688 0 : Complex * wghtSpecStor=wghtSpec.getStorage(weightCopy);
1689 0 : Bool * flagStor=flag.getStorage(flagCopy);
1690 0 : Bool * rowFlagStor=rowFlag.getStorage(rowFlagCopy);
1691 0 : Double * uvwStor=uvw.getStorage(uvwCopy);
1692 0 : Int* ant1Stor=ant1.getStorage(ant1Copy);
1693 0 : Int* ant2Stor=ant2.getStorage(ant2Copy);
1694 0 : Double* timeCenStor= timeCen.getStorage(timeCenCopy);
1695 0 : Double* sumweightStor=sumWeight_p.getStorage(sumWeightCopy);
1696 0 : Double* numvisStor=numVis_p.getStorage(numVisCopy);
1697 : // Vector<Double> invLambda=visFreq/C::c;
1698 :
1699 : //Fill the caches in the master thread
1700 : //// these are thread unsafe...unless already cached
1701 0 : vb.flagRow(); vb.uvw(); vb.antenna1(); vb.antenna2(); vb.time();
1702 0 : vb.nRows(); vb.nCorrelations();
1703 0 : if(hasCorrected)
1704 0 : vb.visCubeCorrected();
1705 : else
1706 0 : vb.visCube();
1707 0 : vb.flagCube();
1708 0 : vb.weight();
1709 : ///////////////////////////
1710 0 : Int nth=1;
1711 : #ifdef _OPENMP
1712 0 : nth=min(nchan_p, omp_get_max_threads());
1713 0 : omp_set_dynamic(0);
1714 : #endif
1715 0 : #pragma omp parallel for firstprivate(refFreq, scale, hasCorrected, needRot, fracbw, gridStor, wghtSpecStor, flagStor, rowFlagStor, uvwStor, ant1Stor, ant2Stor, timeCenStor, sumweightStor, numvisStor ) shared(phasor, visFreq) num_threads(nth) schedule(dynamic, 1)
1716 :
1717 : for(Int outchan=0; outchan < nchan_p; ++outchan){
1718 : //cerr << "outchan " << outchan << " " << chanMapRev_p[outchan] << endl;
1719 :
1720 : multiThrLoop(outchan, vb, refFreq, scale, hasCorrected, needRot, phasor, visFreq,
1721 : fracbw, gridStor, wghtSpecStor, flagStor, rowFlagStor, uvwStor, ant1Stor,
1722 : ant2Stor, timeCenStor, sumweightStor, numvisStor, startchan, endchan, convFunc, convSupport, wScale,
1723 : convSampling);
1724 :
1725 : /*for(uInt nel=0; nel < chanMapRev_p[outchan].nelements(); ++nel ){
1726 : Int chan=chanMapRev_p[outchan][nel];
1727 : for (Int k=0; k < vb.nRows(); ++k){
1728 : if(!vb.flagRow()[k]){
1729 : Int locu, locv, locw,supp, offu, offv;
1730 : {
1731 : locv=Int(Double(ny_p)/2.0+vb.uvw()(1,k)*refFreq*scale(1)+0.5);
1732 : locu=Int(Double(nx_p)/2.0+vb.uvw()(0,k)*refFreq*scale(0)+0.5);
1733 : offv=Int ((Double(locv)-(Double(ny_p)/2.0+vb.uvw()(1,k)*refFreq*scale(1)))*Double(convSampling)+0.5);
1734 : offu=Int ((Double(locu)-(Double(nx_p)/2.0+vb.uvw()(0,k)*refFreq*scale(0)))*Double(convSampling)+0.5);
1735 : locw=Int(sqrt(fabs(wScale*vb.uvw()(2,k)*refFreq/C::c))+0.5);
1736 : supp=locw < convSupport.shape()[0] ? convSupport(locw) :conSupport(convSupport.nelements()-1) ;
1737 : }
1738 :
1739 : Vector<Int> newrow((2*supp+1)*(2*supp+1), -1);
1740 : if(locv < ny_p && locu < nx_p){
1741 : for (Int yy=0; yy< Int(2*supp+1); ++yy){
1742 : Int newlocv=locv+ yy-supp;
1743 : if(newlocv >0 && newlocv < ny_p){
1744 : for (Int xx=0; xx< Int(2*supp+1); ++xx){
1745 : Int newlocu=locu+xx-supp;
1746 : if(newlocu < nx_p && newlocu >0){
1747 : newrow(yy*(2*supp+1)+xx)=newlocv*nx_p+newlocu;
1748 : }
1749 : }
1750 : }
1751 : }
1752 : }
1753 : Complex elphas(1.0, 0.0);
1754 :
1755 : if(chanMap_p(chan) >=startchan && chanMap_p(chan) <=endchan){
1756 : if(needRot){
1757 : Double phasmult=phasor(k)*invLambda(chan);
1758 : Double s, c;
1759 : SINCOS(phasmult, s, c);
1760 : elphas=Complex(c, s);
1761 : }
1762 : if(fracbw > 0.05)
1763 : {
1764 : locv=Int(Double(ny_p)/2.0+vb.uvw()(1,k)*visFreq(chan)*scale(1)+0.5);
1765 : locu=Int(Double(nx_p)/2.0+vb.uvw()(0,k)*visFreq(chan)*scale(0)+0.5);
1766 : }
1767 : if(locv < ny_p && locu < nx_p){
1768 : for (Int yy=0; yy< Int(2*supp+1); ++yy){
1769 : Int locy=abs((yy-supp)*convSampling+offv);
1770 :
1771 : for (Int xx=0; xx< Int(2*supp+1); ++xx){
1772 : Int locx=abs((xx-supp)*convSampling+offu);
1773 :
1774 : Int jj=yy*(2*supp+1)+xx;
1775 : if(newrow[jj] >-1){
1776 : Complex cwt=convFunc(locx, locy, locw);
1777 : if(vb.uvw()(2,k) > 0.0)
1778 : cwt=conj(cwt);
1779 : if(rowFlag(newrow[jj]) && !(vb.flagRow()(k))){
1780 : rowFlag(newrow[jj])=false;
1781 : uvw(2,newrow[jj])=0;
1782 : ant1(newrow[jj])=vb.antenna1()(k);
1783 : ant2(newrow[jj])=vb.antenna2()(k);
1784 : timeCen(newrow[jj])=vb.time()(k);
1785 :
1786 : }
1787 : Int lechan=chanMap_p(chan)-startchan;
1788 : for(Int pol=0; pol < vb.nCorrelations(); ++pol){
1789 : if((!vb.flagCube()(pol,chan, k)) && (polMap_p(pol)>=0) && (vb.weight()(pol,k)>0.0) && (fabs(cwt) > 0.0)){
1790 : // Double newU= Double(locu-nx_p/2)/refFreq/scale(0);
1791 : // Double newV= Double(locv-ny_p/2)/refFreq/scale(1);
1792 : // Double phaseCorr=((newU/vb.uvw()(0,k)-1)+(newV/vb.uvw()(1,k)-1))*vb.uvw()(2,k)*2.0*C::pi*refFreq/C::c;
1793 : Complex toB=hasCorrected ? vb.visCubeCorrected()(pol,chan,k)*vb.weight()(pol,k)*cwt:
1794 : vb.visCube()(pol,chan,k)*vb.weight()(pol,k)*cwt;
1795 : if(needRot)
1796 : toB *=elphas;
1797 : // Double s, c;
1798 : //SINCOS(phaseCorr, s, c);
1799 : // toB=toB*Complex(c,s);
1800 : //Float elwgt=vb.weight()(pol,k)* fabs(real(cwt));
1801 : //////////////TESTING
1802 : Float elwgt=vb.weight()(pol,k)* real(cwt);
1803 : ///////////////////////////
1804 : grid(polMap_p(pol),lechan, newrow[jj])
1805 : = (grid(polMap_p(pol),lechan, newrow[jj])// *wghtSpec(polMap_p(pol),lechan,newrow[jj])
1806 : + toB); ///(elwgt+wghtSpec(polMap_p(pol),lechan,newrow[jj]));
1807 : flag(polMap_p(pol), lechan, newrow[jj])=false;
1808 : //cerr << "weights " << max(vb.weight()) << " spec " << max(vb.weightSpectrum()) << endl;
1809 : //wghtSpec(polMap_p(pol),chanMap_p(chan), newrow)+=vb.weightSpectrum()(pol, chan, k);
1810 : wghtSpec(polMap_p(pol),lechan, newrow[jj]) += elwgt;
1811 : }
1812 : ///We should do that at the end totally
1813 : //wght(pol,newrow)=median(wghtSpec.xyPlane(newrow).row(pol));
1814 : }
1815 : }//if(newrow)
1816 : }
1817 : }
1818 : }//locu && locv
1819 : }
1820 : }
1821 : }
1822 : //sum wgtspec along channels for weight
1823 : }
1824 : */
1825 : }
1826 0 : grid.putStorage(gridStor, gridCopy);
1827 0 : wghtSpec.putStorage(wghtSpecStor, weightCopy);
1828 0 : flag.putStorage(flagStor, flagCopy);
1829 0 : rowFlag.putStorage(rowFlagStor, rowFlagCopy);
1830 0 : uvw.putStorage(uvwStor,uvwCopy);
1831 0 : ant1.putStorage(ant1Stor,ant1Copy);
1832 0 : ant2.putStorage(ant2Stor, ant2Copy);
1833 0 : timeCen.putStorage(timeCenStor, timeCenCopy);
1834 0 : sumWeight_p.putStorage(sumweightStor, sumWeightCopy);
1835 0 : numVis_p.putStorage(numvisStor, numVisCopy);
1836 :
1837 :
1838 0 : }
1839 :
1840 0 : void MSUVBin::multiThrLoop(const Int outchan, const vi::VisBuffer2& vb, Double refFreq, Vector<Float> scale, Bool hasCorrected, Bool needRot, const Vector<Double>& phasor, const Vector<Double>& visFreq, const Double& fracbw, Complex*& grid,
1841 : Complex*& wghtSpec,
1842 : Bool*& flag, Bool*& rowFlag, Double*& uvw,
1843 : Int*& ant1, Int*& ant2, Double*& timeCen, Double*& sumWeight, Double*& numVis,
1844 : const Int startchan, const Int endchan,
1845 : const Cube<Complex>& convFunc, const Vector<Int>& convSupport, const Double wScale, const Int convSampling ){
1846 :
1847 0 : for(uInt nel=0; nel < chanMapRev_p[outchan].nelements(); ++nel ){
1848 0 : Int chan=chanMapRev_p[outchan][nel];
1849 0 : Double invLambda=visFreq[chan]/C::c;
1850 0 : for (rownr_t k=0; k < vb.nRows(); ++k){
1851 0 : if(!vb.flagRow()[k]){
1852 : Int locu, locv, locw,supp, offu, offv;
1853 : {
1854 0 : Double centx=Double(nx_p)/2.0;
1855 0 : Double centy=Double(ny_p)/2.0;
1856 0 : locv=Int(centy+vb.uvw()(1,k)*refFreq*scale(1)+0.5);
1857 0 : locu=Int(centx+vb.uvw()(0,k)*refFreq*scale(0)+0.5);
1858 0 : offv=Int ((Double(locv)-(centy+vb.uvw()(1,k)*refFreq*scale(1)))*Double(convSampling)+0.5);
1859 0 : offu=Int ((Double(locu)-(centx+vb.uvw()(0,k)*refFreq*scale(0)))*Double(convSampling)+0.5);
1860 0 : locw=Int(sqrt(fabs(wScale*vb.uvw()(2,k)*refFreq/C::c))+0.5);
1861 0 : supp=locw < convSupport.shape()[0] ? convSupport(locw) :convSupport(convSupport.nelements()-1) ;
1862 : }
1863 :
1864 0 : Vector<Int> newrow((2*supp+1)*(2*supp+1), -1);
1865 0 : if(locv < ny_p && locu < nx_p){
1866 0 : for (Int yy=0; yy< Int(2*supp+1); ++yy){
1867 0 : Int newlocv=locv+ yy-supp;
1868 0 : if(newlocv >0 && newlocv < ny_p){
1869 0 : for (Int xx=0; xx< Int(2*supp+1); ++xx){
1870 0 : Int newlocu=locu+xx-supp;
1871 0 : if(newlocu < nx_p && newlocu >0){
1872 0 : newrow(yy*(2*supp+1)+xx)=newlocv*nx_p+newlocu;
1873 : }
1874 : }
1875 : }
1876 : }
1877 : }
1878 0 : Complex elphas(1.0, 0.0);
1879 :
1880 0 : if(chanMap_p(chan) >=startchan && chanMap_p(chan) <=endchan){
1881 0 : if(needRot){
1882 0 : Double phasmult=phasor(k)*invLambda;
1883 : Double s, c;
1884 0 : SINCOS(phasmult, s, c);
1885 0 : elphas=Complex(c, s);
1886 : }
1887 0 : if(fracbw > 0.05)
1888 : {
1889 0 : locv=Int(Double(ny_p)/2.0+vb.uvw()(1,k)*visFreq(chan)*scale(1)+0.5);
1890 0 : locu=Int(Double(nx_p)/2.0+vb.uvw()(0,k)*visFreq(chan)*scale(0)+0.5);
1891 : }
1892 0 : if(locv < ny_p && locu < nx_p){
1893 0 : for (Int yy=0; yy< Int(2*supp+1); ++yy){
1894 0 : Int locy=abs((yy-supp)*convSampling+offv);
1895 :
1896 0 : for (Int xx=0; xx< Int(2*supp+1); ++xx){
1897 0 : Int locx=abs((xx-supp)*convSampling+offu);
1898 :
1899 0 : Int jj=yy*(2*supp+1)+xx;
1900 : ////////TEST
1901 : //Int CENTER=2*supp*(supp+1);
1902 : //jj=CENTER;
1903 : /////////////////End TEST
1904 0 : if(newrow[jj] >-1){
1905 0 : Complex cwt=convFunc(locx, locy, locw);
1906 : //Float fracconv=fabs(cwt)/fabs(convFunc(0,0,locw));
1907 : /////TEST
1908 : //cwt=1.0;
1909 : //////
1910 0 : if(vb.uvw()(2,k) > 0.0)
1911 0 : cwt=conj(cwt);
1912 : //if(rowFlag[newrow[jj]] && !(vb.flagRow()(k))){
1913 0 : if(!(vb.flagRow()(k))){
1914 0 : rowFlag[newrow[jj]]=false;
1915 0 : if(yy==supp and xx==supp){
1916 0 : uvw[2+newrow[jj]*3]=0;
1917 0 : for(Int pol=0; pol < vb.nCorrelations(); ++pol){
1918 0 : numVis[npol_p*chanMap_p(chan)+polMap_p(pol)] += 1.0;
1919 0 : sumWeight[npol_p*chanMap_p(chan)+polMap_p(pol)] += vb.weight()(pol,k) ;
1920 : }
1921 : }
1922 : /*else{
1923 : if(uvw[2+newrow[jj]*3] !=0)
1924 : uvw[2+newrow[jj]*3]=-666;
1925 : }*/
1926 0 : ant1[newrow[jj]]=vb.antenna1()(k);
1927 0 : ant2[newrow[jj]]=vb.antenna2()(k);
1928 0 : timeCen[newrow[jj]]=vb.time()(k);
1929 :
1930 : }
1931 0 : Int lechan=chanMap_p(chan)-startchan;
1932 0 : for(Int pol=0; pol < vb.nCorrelations(); ++pol){
1933 0 : if((!vb.flagCube()(pol,chan, k)) && (polMap_p(pol)>=0) /*&& (vb.weight()(pol,k)>0.0) && fabs(cwt) > 0.0//// (fracconv > 5e-2)*/){
1934 : // Double newU= Double(locu-nx_p/2)/refFreq/scale(0);
1935 : // Double newV= Double(locv-ny_p/2)/refFreq/scale(1);
1936 : // Double phaseCorr=((newU/vb.uvw()(0,k)-1)+(newV/vb.uvw()(1,k)-1))*vb.uvw()(2,k)*2.0*C::pi*refFreq/C::c;
1937 : Complex toB=hasCorrected ?
1938 0 : vb.visCubeCorrected()(pol,chan,k)*vb.weight()(pol,k)*cwt :
1939 0 : vb.visCube()(pol,chan,k)*vb.weight()(pol,k)*cwt;
1940 0 : if(needRot)
1941 0 : toB *=elphas;
1942 : // Double s, c;
1943 : //SINCOS(phaseCorr, s, c);
1944 : // toB=toB*Complex(c,s);
1945 : //Float elwgt=vb.weight()(pol,k)* fabs(real(cwt));
1946 : //////////////TESTING
1947 0 : Complex elwgt=vb.weight()(pol,k)* (cwt);
1948 : ///////////////////////////
1949 0 : ooLong cubindx=ooLong(newrow[jj])*ooLong(npol_p)*ooLong(endchan-startchan+1)+ooLong(lechan*npol_p)+ooLong(polMap_p(pol));
1950 : /* if(newrow[jj]==2023010 && outchan==0){
1951 : cerr << jj << " newrow[jj] " << newrow[jj] << " polMap_p(pol) " <<polMap_p(pol) << " cubindex " << cubindx << endl;
1952 : cerr << "uvw " << vb.uvw().column(k) << " weight " << vb.weight()(pol,k) << " elwgt " << elwgt << " cwt " << cwt << endl;
1953 : cerr << "revmap " << chanMapRev_p[outchan].nelements() << " support " << supp << " locu locv " << locu << " " <<locv << endl;
1954 : }*/
1955 0 : grid[cubindx]
1956 0 : = (grid[cubindx]// *wghtSpec(polMap_p(pol),lechan,newrow[jj])
1957 0 : + toB); ///(elwgt+wghtSpec(polMap_p(pol),lechan,newrow[jj]));
1958 0 : flag[cubindx]=false;
1959 : //cerr << "weights " << max(vb.weight()) << " spec " << max(vb.weightSpectrum()) << endl;
1960 : //wghtSpec(polMap_p(pol),chanMap_p(chan), newrow)+=vb.weightSpectrum()(pol, chan, k);
1961 0 : wghtSpec[cubindx] += elwgt;
1962 : //wghtSpec[cubindx] +=vb.weight()(pol,k)/(4*supp*supp);
1963 : }
1964 : ///We should do that at the end totally
1965 : //wght(pol,newrow)=median(wghtSpec.xyPlane(newrow).row(pol));
1966 : }
1967 : }//if(newrow)
1968 : }
1969 : }
1970 : }//locu && locv
1971 : }
1972 0 : }
1973 : }
1974 : //sum wgtspec along channels for weight
1975 : /*for (Int pol=0; pol < wght.shape()(0); ++pol){
1976 : //cerr << "shape min max "<< median(wghtSpec.xyPlane(newrow).row(pol)) << " " << min(wghtSpec.xyPlane(newrow).row(pol)) << " "<< max(wghtSpec.xyPlane(newrow).row(pol)) << endl;
1977 : wght(pol,newrow)=median(wghtSpec.xyPlane(newrow).row(pol));
1978 : //cerr << "pol " << pol << " newrow "<< newrow << " weight "<< wght(pol, newrow) << endl;
1979 : }
1980 : */
1981 : }
1982 :
1983 0 : }
1984 :
1985 :
1986 0 : void MSUVBin::locateFlagFromGrid(vi::VisBuffer2& vb, Cube<Bool>& datFlag,
1987 : Cube<Float>& wghtSpec,
1988 : Cube<Bool>& flag, Vector<Bool>& rowFlag, Matrix<Double>& uvw, Vector<Int>& ant1,
1989 : Vector<Int>& ant2, Vector<Double>& timeCen, const Int startchan, const Int endchan){
1990 : //all pixel that are touched the flag and flag Row shall be unset and the w be assigned
1991 : //later we'll deal with multiple w for the same uv
1992 : //we need polmap and chanmap;
1993 :
1994 : Double fracbw;
1995 0 : if(!datadescMap(vb, fracbw)) return;
1996 : //cerr << "fracbw " << fracbw << endl;
1997 0 : SpectralCoordinate spec=csys_p.spectralCoordinate(2);
1998 0 : DirectionCoordinate thedir=csys_p.directionCoordinate(0);
1999 0 : Double refFreq=SpectralImageUtil::worldFreq(csys_p, Double(nchan_p/2));
2000 : //Double refFreq=SpectralImageUtil::worldFreq(csys_p, Double(0));
2001 0 : Vector<Float> scale(2);
2002 0 : scale(0)=fabs(nx_p*thedir.increment()(0))/C::c;
2003 0 : scale(1)=fabs(ny_p*thedir.increment()(1))/C::c;
2004 : //Dang i thought the new vb will return Data or FloatData if correctedData was
2005 : //not there
2006 : //Bool hasCorrected=!(MSMainColumns(vb.ms()).correctedData().isNull());
2007 : //locateuvw(locuv, vb.uvw());
2008 0 : Vector<Double> visFreq=vb.getFrequencies(0, MFrequency::LSRK);
2009 0 : for (rownr_t k=0; k < vb.nRows(); ++k){
2010 0 : if(!vb.flagRow()[k]){
2011 : Int locu, locv;
2012 : {
2013 0 : locv=Int(Double(ny_p)/2.0+vb.uvw()(1,k)*refFreq*scale(1)+0.5);
2014 0 : locu=Int(Double(nx_p)/2.0+vb.uvw()(0,k)*refFreq*scale(0)+0.5);
2015 :
2016 : }
2017 :
2018 0 : for(Int chan=0; chan < vb.nChannels(); ++chan ){
2019 0 : if(chanMap_p(chan) >=startchan && chanMap_p(chan) <=endchan){
2020 : //Double outChanFreq;
2021 : //spec.toWorld(outChanFreq, Double(chanMap_p(chan)));
2022 0 : if(fracbw > 0.05)
2023 : {
2024 0 : locv=Int(Double(ny_p)/2.0+vb.uvw()(1,k)*visFreq(chan)*scale(1)+0.5);
2025 0 : locu=Int(Double(nx_p)/2.0+vb.uvw()(0,k)*visFreq(chan)*scale(0)+0.5);
2026 : }
2027 0 : if(locv < ny_p && locu < nx_p && locv >=0 && locu >=0){
2028 0 : Int newrow=locv*nx_p+locu;
2029 0 : if(rowFlag(newrow) && !(vb.flagRow()(k))){
2030 0 : rowFlag(newrow)=false;
2031 : /////TEST
2032 : //uvw(0,newrow)=vb.uvw()(0,k);
2033 : //uvw(1,newrow)=vb.uvw()(1,k);
2034 : /////
2035 0 : uvw(2,newrow)=vb.uvw()(2,k);
2036 : //cerr << newrow << " rowids " << vb.rowIds()[k] << " uvw2 " << uvw(2 ,newrow) << endl;
2037 0 : ant1(newrow)=vb.antenna1()(k);
2038 0 : ant2(newrow)=vb.antenna2()(k);
2039 0 : timeCen(newrow)=vb.time()(k);
2040 :
2041 : }
2042 0 : for(Int pol=0; pol < vb.nCorrelations(); ++pol){
2043 0 : if((!vb.flagCube()(pol,chan, k)) && (polMap_p(pol)>=0) && (vb.weight()(pol,k)>0.0)){
2044 : // Double newU= Double(locu-nx_p/2)/refFreq/scale(0);
2045 : // Double newV= Double(locv-ny_p/2)/refFreq/scale(1);
2046 : // Double phaseCorr=((newU/vb.uvw()(0,k)-1)+(newV/vb.uvw()(1,k)-1))*vb.uvw()(2,k)*2.0*C::pi*refFreq/C::c;
2047 : // Complex toB=hasCorrected ? vb.visCubeCorrected()(pol,chan,k)*vb.weight()(pol,k):
2048 : // vb.visCube()(pol,chan,k)*vb.weight()(pol,k);
2049 : // Double s, c;
2050 : //SINCOS(phaseCorr, s, c);
2051 : // toB=toB*Complex(c,s);
2052 : // grid(polMap_p(pol),chanMap_p(chan)-startchan, newrow)
2053 : // = (grid(polMap_p(pol),chanMap_p(chan)-startchan, newrow)
2054 : // + toB); ///(vb.weight()(pol,k)+wghtSpec(polMap_p(pol),chanMap_p(chan)-startchan,newrow));
2055 0 : if(flag(polMap_p(pol),chanMap_p(chan)-startchan, newrow) || wghtSpec(polMap_p(pol),chanMap_p(chan)-startchan, newrow)==0.0) {
2056 0 : datFlag(pol,chan,k)=true;
2057 :
2058 : }
2059 : //cerr << "weights " << max(vb.weight()) << " spec " << max(vb.weightSpectrum()) << endl;
2060 : //wghtSpec(polMap_p(pol),chanMap_p(chan), newrow)+=vb.weightSpectrum()(pol, chan, k);
2061 : //wghtSpec(polMap_p(pol),chanMap_p(chan)-startchan, newrow) += vb.weight()(pol,k);
2062 : }
2063 : ///We should do that at the end totally
2064 : //wght(pol,newrow)=median(wghtSpec.xyPlane(newrow).row(pol));
2065 : }
2066 : }//locu && locv
2067 : }
2068 : }
2069 : //sum wgtspec along channels for weight
2070 : /*for (Int pol=0; pol < wght.shape()(0); ++pol){
2071 : //cerr << "shape min max "<< median(wghtSpec.xyPlane(newrow).row(pol)) << " " << min(wghtSpec.xyPlane(newrow).row(pol)) << " "<< max(wghtSpec.xyPlane(newrow).row(pol)) << endl;
2072 : wght(pol,newrow)=median(wghtSpec.xyPlane(newrow).row(pol));
2073 : //cerr << "pol " << pol << " newrow "<< newrow << " weight "<< wght(pol, newrow) << endl;
2074 : }
2075 : */
2076 : }
2077 : }
2078 :
2079 :
2080 :
2081 :
2082 0 : }
2083 0 : Bool MSUVBin::saveData(const Cube<Complex>& grid, const Cube<Bool>&flag, const Vector<Bool>& rowFlag,
2084 : const Cube<Float>&wghtSpec,
2085 : const Matrix<Double>& uvw, const Vector<Int>& ant1,
2086 : const Vector<Int>& ant2, const Vector<Double>& timeCen, const Int startchan, const Int endchan, const Cube<Float>& imagWghtSpec){
2087 0 : Bool retval=true;
2088 0 : MSColumns msc(*outMsPtr_p);
2089 0 : if(!existOut_p && startchan==0){
2090 0 : fillSubTables();
2091 : // msc.uvw().putColumn(uvw);
2092 :
2093 : }
2094 0 : msc.uvw().putColumn(uvw);
2095 0 : Slicer elslice(IPosition(2, 0, startchan), IPosition(2,npol_p, endchan-startchan+1));
2096 : //lets put ny_p slices
2097 0 : Int nchan=endchan-startchan+1;
2098 0 : Slice polslice(0,npol_p);
2099 0 : Slice chanslice(0,nchan);
2100 0 : for (Int k=0; k <ny_p; ++k){
2101 : //Slicer rowslice(IPosition(1, k*nx_p), IPosition(1,nx_p));
2102 0 : RefRows rowslice(k*nx_p, (k+1)*nx_p-1);
2103 0 : msc.data().putColumnCells(rowslice, elslice, grid(polslice, chanslice, Slice(k*nx_p, nx_p)));
2104 : //msc.data().putColumnRange(rowslice, elslice, grid(polslice, chanslice, Slice(k*nx_p, nx_p)));
2105 0 : msc.weightSpectrum().putColumnCells(rowslice, elslice, wghtSpec(polslice,chanslice, Slice(k*nx_p, nx_p)));
2106 : //msc.weightSpectrum().putColumnRange(rowslice, elslice, wghtSpec(polslice,chanslice, Slice(k*nx_p, nx_p)));
2107 : //////msc.weight().putColumn(wght);
2108 0 : msc.flag().putColumnCells(rowslice, elslice,flag(polslice, chanslice, Slice(k*nx_p,nx_p)));
2109 : //msc.flag().putColumnRange(rowslice, elslice,flag(polslice, chanslice, Slice(k*nx_p,nx_p)));
2110 0 : }
2111 0 : if(doW_p){
2112 0 : ArrayColumn<Float> imagWgtSpecCol(*outMsPtr_p, "IMAGINARY_WEIGHT_SPECTRUM");
2113 0 : for (Int k=0; k <ny_p; ++k){
2114 0 : RefRows rowslice(k*nx_p, (k+1)*nx_p-1);
2115 0 : imagWgtSpecCol.putColumnCells(rowslice, elslice, imagWghtSpec(polslice,chanslice, Slice(k*nx_p, nx_p)));
2116 0 : }
2117 0 : }
2118 0 : if(endchan==nchan_p-1){
2119 0 : msc.flagRow().putColumn(rowFlag);
2120 0 : msc.antenna1().putColumn(ant1);
2121 0 : msc.antenna2().putColumn(ant2);
2122 0 : Cube<Float> spectralweight;
2123 0 : msc.weightSpectrum().getColumn(spectralweight);
2124 0 : Matrix<Float> weight(npol_p, wghtSpec.shape()[2]);
2125 0 : for (Int row=0; row < wghtSpec.shape()[2]; ++row){
2126 0 : for (Int pol=0; pol < npol_p; ++pol){
2127 : //cerr << "shape min max "<< median(wghtSpec.xyPlane(newrow).row(pol)) << " " << min(wghtSpec.xyPlane(newrow).row(pol)) << " "<< max(wghtSpec.xyPlane(newrow).row(pol)) << endl;
2128 0 : weight(pol,row)=min(spectralweight.xyPlane(row).row(pol));
2129 : //weight(pol,row)=1.0;
2130 : //if(!rowFlag(row))
2131 : // cerr << "pol " << pol << " row "<< row << " median "<< weight(pol, row) << " min-max " << (min(wghtSpec.xyPlane(row).row(pol))+max(wghtSpec.xyPlane(row).row(pol)))/2.0 << " mean " << mean(wghtSpec.xyPlane(row).row(pol)) << endl;
2132 : }
2133 : }
2134 :
2135 0 : msc.weight().putColumn(weight);
2136 0 : Matrix<Float> sigma(weight.shape());
2137 0 : for (uInt k=0; k<weight.shape()[1]; ++k)
2138 0 : for (uInt j=0; j<weight.shape()[0]; ++j)
2139 0 : sigma(j,k)=weight(j,k)>0.0 ? 1/sqrt(weight(j,k)) :0.0;
2140 0 : msc.sigma().putColumn(sigma);
2141 :
2142 0 : msc.time().putColumn(timeCen);
2143 0 : msc.timeCentroid().putColumn(timeCen);
2144 0 : }
2145 0 : return retval;
2146 :
2147 :
2148 0 : }
2149 0 : Bool MSUVBin::saveData(const Cube<Complex>& grid, const Cube<Bool>&flag, const Vector<Bool>& rowFlag,
2150 : const Cube<Float>&wghtSpec, const Matrix<Float>& /*wght*/,
2151 : const Matrix<Double>& uvw, const Vector<Int>& ant1,
2152 : const Vector<Int>& ant2, const Vector<Double>& timeCen){
2153 0 : Bool retval=true;
2154 0 : MSColumns msc(*outMsPtr_p);
2155 0 : if(!existOut_p){
2156 0 : fillSubTables();
2157 0 : msc.uvw().putColumn(uvw);
2158 :
2159 : }
2160 :
2161 :
2162 0 : msc.data().putColumn(grid);
2163 0 : msc.weightSpectrum().putColumn(wghtSpec);
2164 : //msc.weight().putColumn(wght);
2165 0 : msc.flag().putColumn(flag);
2166 0 : msc.flagRow().putColumn(rowFlag);
2167 0 : msc.antenna1().putColumn(ant1);
2168 0 : msc.antenna2().putColumn(ant2);
2169 0 : Matrix<Float> weight(npol_p, wghtSpec.shape()[2]);
2170 0 : for (Int row=0; row < wghtSpec.shape()[2]; ++row){
2171 0 : for (Int pol=0; pol < npol_p; ++pol){
2172 : //cerr << "shape min max "<< median(wghtSpec.xyPlane(newrow).row(pol)) << " " << min(wghtSpec.xyPlane(newrow).row(pol)) << " "<< max(wghtSpec.xyPlane(newrow).row(pol)) << endl;
2173 0 : weight(pol,row)=max(wghtSpec.xyPlane(row).row(pol));
2174 : //if(!rowFlag(row))
2175 : // cerr << "pol " << pol << " row "<< row << " median "<< weight(pol, row) << " min-max " << (min(wghtSpec.xyPlane(row).row(pol))+max(wghtSpec.xyPlane(row).row(pol)))/2.0 << " mean " << mean(wghtSpec.xyPlane(row).row(pol)) << endl;
2176 : }
2177 : }
2178 0 : msc.weight().putColumn(weight);
2179 0 : Matrix<Float> sigma(weight.shape());
2180 0 : for (uInt k=0; k<weight.shape()[1]; ++k)
2181 0 : for (uInt j=0; j<weight.shape()[0]; ++j)
2182 0 : sigma(j,k)=weight(j,k)>0.0 ? 1/sqrt(weight(j,k)) :0.0;
2183 0 : msc.sigma().putColumn(sigma);
2184 0 : msc.time().putColumn(timeCen);
2185 0 : msc.timeCentroid().putColumn(timeCen);
2186 :
2187 0 : return retval;
2188 :
2189 :
2190 0 : }
2191 0 : void MSUVBin::fillSubTables(){
2192 0 : fillFieldTable();
2193 0 : copySubtable("SPECTRAL_WINDOW", mss_p[0]->spectralWindow(), true);
2194 0 : copySubtable("POLARIZATION", mss_p[0]->polarization(), true);
2195 0 : copySubtable("DATA_DESCRIPTION", mss_p[0]->dataDescription(), true);
2196 0 : copySubtable("FEED", mss_p[0]->feed(), false);
2197 0 : copySubtable("OBSERVATION", mss_p[0]->observation(), false);
2198 0 : copySubtable("ANTENNA", mss_p[0]->antenna(), false);
2199 0 : fillDDTables();
2200 :
2201 0 : }
2202 0 : void MSUVBin::fillDDTables(){
2203 : {
2204 0 : MSPolarizationColumns mspol(outMsPtr_p->polarization());
2205 :
2206 0 : if(outMsPtr_p->polarization().nrow()==0)
2207 0 : outMsPtr_p->polarization().addRow();
2208 0 : mspol.numCorr().put(0, npol_p);
2209 0 : mspol.corrType().put(0,whichStokes_p);
2210 0 : Matrix<Int> corrProd(2,npol_p);
2211 0 : corrProd.set(0);
2212 0 : if(npol_p==2){
2213 0 : corrProd(0,1)=1;
2214 0 : corrProd(1,1)=1;
2215 : }
2216 0 : if(npol_p==4){
2217 0 : corrProd(0,1)=0;
2218 0 : corrProd(1,1)=1;
2219 0 : corrProd(0,2)=1;
2220 0 : corrProd(1,2)=0;
2221 0 : corrProd(0,3)=1;
2222 0 : corrProd(1,3)=1;
2223 : }
2224 0 : mspol.corrProduct().put(0,corrProd);
2225 0 : mspol.flagRow().put(0,false);
2226 0 : }
2227 : ///Now with Spectral window
2228 : {
2229 0 : MSSpWindowColumns msSpW(outMsPtr_p->spectralWindow());
2230 0 : if(outMsPtr_p->spectralWindow().nrow()==0){
2231 0 : MSTransformDataHandler::addOptionalColumns(mss_p[0]->spectralWindow(),
2232 0 : outMsPtr_p->spectralWindow());
2233 0 : outMsPtr_p->spectralWindow().addRow();
2234 : }
2235 0 : msSpW.name().put(0,"none");
2236 0 : msSpW.ifConvChain().put(0,0);
2237 0 : msSpW.numChan().put(0,nchan_p);
2238 0 : SpectralCoordinate spec=csys_p.spectralCoordinate(2);
2239 : Double refFreq;
2240 0 : spec.toWorld(refFreq,0.0);
2241 0 : Double chanBandWidth=spec.increment()[0];
2242 0 : Vector<Double> chanFreq(nchan_p),resolution(nchan_p);
2243 0 : for (Int i=0; i < nchan_p; ++i) {
2244 0 : spec.toWorld(chanFreq[i], Double(i));
2245 : }
2246 0 : resolution=chanBandWidth;
2247 0 : msSpW.chanFreq().put(0,chanFreq);
2248 0 : msSpW.chanWidth().put(0,resolution);
2249 0 : msSpW.effectiveBW().put(0,resolution);
2250 0 : msSpW.refFrequency().put(0,refFreq);
2251 0 : msSpW.resolution().put(0,resolution);
2252 0 : msSpW.totalBandwidth().put(0,abs(nchan_p*chanBandWidth));
2253 0 : msSpW.netSideband().put(0,1);
2254 0 : msSpW.freqGroup().put(0,0);
2255 0 : msSpW.freqGroupName().put(0,"none");
2256 0 : msSpW.flagRow().put(0,false);
2257 0 : msSpW.measFreqRef().put(0,MFrequency::LSRK);
2258 0 : }
2259 : //Now the DD
2260 : {
2261 0 : MSDataDescColumns msDD(outMsPtr_p->dataDescription());
2262 0 : if(outMsPtr_p->dataDescription().nrow()==0)
2263 0 : outMsPtr_p->dataDescription().addRow();
2264 0 : msDD.spectralWindowId().put(0,0);
2265 0 : msDD.polarizationId().put(0,0);
2266 0 : msDD.flagRow().put(0,false);
2267 0 : }
2268 0 : MSColumns(*outMsPtr_p).dataDescId().fillColumn(0);
2269 :
2270 :
2271 0 : }
2272 0 : void MSUVBin::copySubtable(const String& tabName, const Table& inTab,const Bool norows)
2273 : {
2274 : //outMsPtr_p->closeSubTables();
2275 0 : String outName(outMsPtr_p->tableName() + '/' + tabName);
2276 :
2277 0 : if (PlainTable::tableCache()(outName)){
2278 : //cerr << "cpy subtable "<< outName << endl;
2279 0 : Table outTab(outName, Table::Update);
2280 0 : if(norows){
2281 0 : Vector<rownr_t> rownums=outTab.rowNumbers();
2282 0 : outTab.removeRow(rownums);
2283 0 : }
2284 : else{
2285 0 : TableCopy::copySubTables(outTab, inTab);
2286 0 : TableCopy::copyInfo(outTab, inTab);
2287 : //cerr << "ROWS "<< inTab.nrow() << " " << outTab.nrow() << endl;
2288 0 : TableCopy::copyRows(outTab, inTab);
2289 : /*outTab=Table();
2290 : cerr << "copying " << tabName << endl;
2291 : inTab.copy(outName, Table::New);
2292 : */
2293 : }
2294 0 : }
2295 : else{
2296 0 : inTab.deepCopy(outName, Table::New, false, Table::AipsrcEndian, norows);
2297 : }
2298 0 : Table outTab(outName, Table::Update);
2299 0 : outMsPtr_p->rwKeywordSet().defineTable(tabName, outTab);
2300 0 : outMsPtr_p->initRefs();
2301 :
2302 0 : return;
2303 0 : }
2304 :
2305 0 : void MSUVBin::fillFieldTable() {
2306 0 : MSFieldColumns msField(outMsPtr_p->field());
2307 :
2308 0 : if(outMsPtr_p->field().nrow()==0)
2309 0 : outMsPtr_p->field().addRow();
2310 :
2311 0 : msField.sourceId().put(0, 0);
2312 0 : msField.code().put(0, "");
2313 0 : String sourceName="MSUVBIN";
2314 0 : Int fieldId=MSColumns(*mss_p[0]).fieldId()(0);
2315 :
2316 0 : sourceName=MSFieldColumns((mss_p[0])->field()).name()(fieldId);
2317 0 : msField.name().put(0, sourceName);
2318 0 : Int numPoly = 0;
2319 :
2320 : //
2321 0 : Vector<MDirection> radecMeas(1);
2322 0 : radecMeas(0)=phaseCenter_p;
2323 :
2324 0 : Double obsTime=MSColumns(*mss_p[0]).time()(0);
2325 0 : msField.time().put(0, obsTime);
2326 0 : msField.numPoly().put(0, numPoly);
2327 0 : msField.delayDirMeasCol().put(0, radecMeas);
2328 0 : msField.phaseDirMeasCol().put(0, radecMeas);
2329 0 : msField.referenceDirMeasCol().put(0, radecMeas);
2330 0 : msField.flagRow().put(0, false);
2331 0 : MSColumns(*outMsPtr_p).fieldId().fillColumn(0);
2332 0 : }
2333 :
2334 0 : Bool MSUVBin::String2MDirection(const String& theString,
2335 : MDirection& theMeas, const String msname){
2336 :
2337 0 : istringstream istr(theString);
2338 : Int fieldid;
2339 0 : istr >> fieldid;
2340 0 : if(!istr.fail() && msname != ""){
2341 : //We'll interprete string as a field id of ms
2342 0 : MeasurementSet thems(msname);
2343 0 : theMeas=MSFieldColumns(thems.field()).phaseDirMeas(fieldid);
2344 0 : return true;
2345 0 : }
2346 :
2347 :
2348 0 : QuantumHolder qh;
2349 0 : String error;
2350 :
2351 0 : Vector<String> str;
2352 : //In case of compound strings with commas or empty space
2353 0 : sepCommaEmptyToVectorStrings(str, theString);
2354 :
2355 0 : if(str.nelements()==3){
2356 0 : qh.fromString(error, str[1]);
2357 0 : casacore::Quantity val1=qh.asQuantity();
2358 0 : qh.fromString(error, str[2]);
2359 0 : casacore::Quantity val2=qh.asQuantity();
2360 : MDirection::Types tp;
2361 0 : if(!MDirection::getType(tp, str[0])){
2362 0 : ostringstream oss;
2363 0 : oss << "Could not understand Direction frame...defaulting to J2000 " ;
2364 0 : cerr << oss.str() << endl;
2365 0 : tp=MDirection::J2000;
2366 0 : }
2367 0 : theMeas=MDirection(val1,val2, tp);
2368 0 : return true;
2369 0 : }
2370 0 : else if(str.nelements()==2){
2371 0 : qh.fromString(error, str[0]);
2372 0 : casacore::Quantity val1=qh.asQuantity();
2373 0 : qh.fromString(error, str[1]);
2374 0 : casacore::Quantity val2=qh.asQuantity();
2375 0 : theMeas=MDirection(val1, val2);
2376 0 : return true;
2377 0 : }
2378 0 : else if(str.nelements()==1){
2379 : //Must be a string like sun, moon, jupiter
2380 0 : casacore::Quantity val1(0.0, "deg");
2381 0 : casacore::Quantity val2(90.0, "deg");
2382 0 : theMeas=MDirection(val1, val2);
2383 : MDirection::Types ref;
2384 : Int numAll;
2385 : Int numExtra;
2386 : const uInt *dum;
2387 0 : const String *allTypes=MDirection::allMyTypes(numAll, numExtra, dum);
2388 : //if it is SUN moon etc
2389 0 : if(MDirection::getType(ref,str[0])){
2390 :
2391 0 : theMeas=MDirection(val1, val2, ref);
2392 0 : return true;
2393 : }
2394 0 : if(MeasTable::Source(theMeas, str[0])){
2395 0 : return true;
2396 : }
2397 0 : if(!MDirection::getType(ref, str[0])){
2398 0 : Vector<String> all(numExtra);
2399 0 : for(Int k =0; k < numExtra; ++k){
2400 0 : all[k]=*(allTypes+numAll-k-1);
2401 : }
2402 0 : ostringstream oss;
2403 0 : oss << "Could not understand Direction string " <<str[0] << "\n" ;
2404 0 : oss << "Valid ones are " << all;
2405 0 : cerr << oss.str() << " or one of the valid known sources in the data repos" << endl;
2406 0 : theMeas=MDirection(val1, val2);
2407 0 : return false;
2408 0 : }
2409 :
2410 0 : }
2411 :
2412 :
2413 :
2414 :
2415 : ///If i am here i don't know how to interprete this
2416 :
2417 :
2418 0 : return false;
2419 0 : }
2420 :
2421 :
2422 0 : void MSUVBin::makeSFConv(Cube<Complex>& convFunc, Vector<Int>& convSupport, Double& wScale, Int& convSampling, Int& convSize){
2423 0 : Vector<Double> uvOffset(2);
2424 0 : uvOffset(0)=Double(nx_p)/2.0;
2425 0 : uvOffset(1)=Double(ny_p)/2.0;
2426 0 : Vector<Double> uvScale(2);
2427 : //cerr << "increments " << csys_p.increment() << endl;
2428 0 : uvScale(0)=fabs(nx_p*csys_p.increment()(0));
2429 0 : uvScale(1)=fabs(ny_p*csys_p.increment()(1));
2430 : /////no w
2431 0 : wScale=0.0;
2432 :
2433 : ConvolveGridder<Double, Complex>
2434 0 : gridder(IPosition(2, nx_p, ny_p), uvScale, uvOffset, "SF");
2435 0 : convSupport.resize(1);
2436 0 : convSupport=gridder.cSupport()(0);
2437 0 : convSize=convSupport(0);
2438 0 : convSampling=gridder.cSampling();
2439 : //cerr << " support " << gridder.cSupport() << " sampling " << convSampling << endl;
2440 :
2441 0 : Vector<Double> cfunc1D=gridder.cFunction();
2442 :
2443 : //cerr << "convFunc " << cfunc1D << endl;
2444 0 : convFunc.resize(convSampling*convSupport(0), convSampling*convSupport(0), 1);
2445 0 : convFunc.set(Complex(0.0));
2446 0 : for (Int k=0; k < convSampling*convSupport(0); ++k){
2447 0 : for (int j=0; j < convSampling*convSupport(0); ++j){
2448 0 : convFunc(j, k,0)=Complex(cfunc1D(k)*cfunc1D(j));
2449 : }
2450 : }
2451 :
2452 0 : }
2453 :
2454 0 : void MSUVBin::makeWConv(vi::VisibilityIterator2& iter, Cube<Complex>& convFunc, Vector<Int>& convSupport,
2455 : Double& wScale, Int& convSampling, Int& convSize){
2456 : ///////Let's find ... the maximum w
2457 0 : vi::VisBuffer2 *vb=iter.getVisBuffer();
2458 0 : Double maxW=0.0;
2459 0 : Double minW=1e99;
2460 0 : Double nval=0;
2461 0 : Double rmsW=0.0;
2462 0 : for (iter.originChunks(); iter.moreChunks(); iter.nextChunk()) {
2463 0 : for (iter.origin(); iter.more(); iter.next()) {
2464 0 : maxW=max(maxW, max(abs(vb->uvw().row(2)*max(vb->getFrequencies(0))))/C::c);
2465 0 : minW=min(minW, min(abs(vb->uvw().row(2)*min(vb->getFrequencies(0))))/C::c);
2466 : ///////////some shenanigans as some compilers is confusing * operator for vector
2467 0 : Vector<Double> elw;
2468 0 : elw=vb->uvw().row(2);
2469 0 : elw*=vb->uvw().row(2);
2470 : //////////////////////////////////////////////////
2471 0 : rmsW+=sum(elw);
2472 :
2473 0 : nval+=Double(vb->nRows());
2474 0 : }
2475 : }
2476 0 : rmsW=sqrt(rmsW/Double(nval))*max(vb->getFrequencies(0))/C::c;
2477 0 : Double maxUVW=rmsW < 0.5*(minW+maxW) ? 1.05*maxW: (rmsW /(0.5*((minW)+maxW))*1.05*maxW) ;
2478 : //Int wConvSize=Int(maxUVW*fabs(sin(fabs(csys_p.increment()(0))*max(nx_p, ny_p)/2.0)));
2479 : //////////TEST
2480 0 : Int wConvSize=2*Int(maxUVW*fabs(sin(fabs(csys_p.increment()(0))*max(nx_p, ny_p)/2.0)));
2481 : /////////////////
2482 0 : wScale=Double((wConvSize-1)*(wConvSize-1))/maxUVW;
2483 0 : if(wConvSize <2)
2484 0 : ThrowCc("Should not be using wprojection");
2485 0 : Int maxMemoryMB=HostInfo::memoryTotal(true)/1024;
2486 0 : Double maxConvSizeConsidered=sqrt(Double(maxMemoryMB)/8.0*1024.0*1024.0/Double(wConvSize));
2487 0 : CompositeNumber cn(Int(maxConvSizeConsidered/2.0)*2);
2488 : //cerr << "max ConvSize considered " << maxConvSizeConsidered << endl;
2489 0 : convSampling=10;
2490 0 : convSize=max(nx_p, ny_p);
2491 0 : convSize=min(convSize,(Int)cn.nearestEven(Int(maxConvSizeConsidered/2.0)*2));
2492 0 : Int maxConvSize=convSize;
2493 0 : CoordinateSystem coords=csys_p;
2494 0 : DirectionCoordinate dc=coords.directionCoordinate(0);
2495 0 : Vector<Double> sampling;
2496 0 : sampling = dc.increment();
2497 0 : sampling*=Double(convSampling);
2498 0 : sampling[0]*=Double(cn.nextLargerEven(Int(Float(nx_p)-0.5)))/Double(convSize);
2499 0 : sampling[1]*=Double(cn.nextLargerEven(Int(Float(ny_p)-0.5)))/Double(convSize);
2500 0 : dc.setIncrement(sampling);
2501 0 : Vector<Double> unitVec(2);
2502 0 : unitVec=convSize/2;
2503 0 : dc.setReferencePixel(unitVec);
2504 :
2505 : // Set the reference value to that of the image center for sure.
2506 : {
2507 : // dc.setReferenceValue(mTangent_p.getAngle().getValue());
2508 0 : MDirection wcenter;
2509 0 : Vector<Double> pcenter(2);
2510 0 : pcenter(0) = nx_p/2;
2511 0 : pcenter(1) = ny_p/2;
2512 0 : coords.directionCoordinate(0).toWorld( wcenter, pcenter );
2513 0 : dc.setReferenceValue(wcenter.getAngle().getValue());
2514 0 : }
2515 0 : coords.replaceCoordinate(dc, 0);
2516 0 : IPosition pbShape(4, convSize, convSize, 1, 1);
2517 0 : TempImage<Complex> twoDPB(pbShape, coords);
2518 0 : convFunc.resize(); // break any reference
2519 0 : convFunc.resize(convSize/2-1, convSize/2-1, wConvSize);
2520 0 : convFunc.set(0.0);
2521 0 : Bool convFuncStor=false;
2522 0 : Complex *convFuncPtr=convFunc.getStorage(convFuncStor);
2523 0 : IPosition start(4, 0, 0, 0, 0);
2524 0 : IPosition pbSlice(4, convSize, convSize, 1, 1);
2525 0 : Vector<Complex> maxes(wConvSize);
2526 : Bool maxdel;
2527 0 : Complex* maxptr=maxes.getStorage(maxdel);
2528 0 : Int inner=convSize/convSampling;
2529 : //cerr << "inner " << inner << endl;
2530 0 : Matrix<Complex> corr(inner, inner);
2531 :
2532 0 : Vector<Double> uvOffset(2);
2533 0 : uvOffset(0)=Double(nx_p)/2.0;
2534 0 : uvOffset(1)=Double(ny_p)/2.0;
2535 0 : Vector<Double> uvScale(2);
2536 0 : Vector<Complex> correction(inner);
2537 : //cerr << "increments " << csys_p.increment() << endl;
2538 0 : uvScale(0)=fabs(nx_p*csys_p.increment()(0));
2539 0 : uvScale(1)=fabs(ny_p*csys_p.increment()(1));
2540 :
2541 :
2542 : ConvolveGridder<Double, Complex>
2543 :
2544 0 : ggridder(IPosition(2, inner, inner), uvScale, uvOffset, "SF");
2545 : /////////////////////////////////////////////////////////
2546 : /////////////Testing
2547 : /*correction.resize(2*inner);
2548 : uvScale *=2.0;
2549 : uvOffset *=2.0;
2550 : ConvolveGridder<Double, Complex>
2551 : ggridder(IPosition(2, 2*inner, 2*inner), uvScale, uvOffset, "SF");
2552 : for (Int iy=-inner/2;iy<inner/2;iy++) {
2553 :
2554 : ggridder.correctX1D(correction, iy+inner);
2555 : corr.row(iy+inner/2)=correction(Slice(inner/2, inner));
2556 : }
2557 : */
2558 : /////////////////////////////////
2559 : //////////////////////////////////////////
2560 :
2561 0 : for (Int iy=-inner/2;iy<inner/2;iy++) {
2562 :
2563 0 : ggridder.correctX1D(correction, iy+inner/2);
2564 0 : corr.row(iy+inner/2)=correction;
2565 : }
2566 :
2567 :
2568 : Bool cpcor;
2569 0 : Complex *cor=corr.getStorage(cpcor);
2570 0 : Double s1=sampling(1);
2571 0 : Double s0=sampling(0);
2572 : ///////////Por FFTPack
2573 0 : Vector<Float> wsave(2*convSize*convSize+15);
2574 0 : Int lsav=2*convSize*convSize+15;
2575 : Bool wsavesave;
2576 0 : Float *wsaveptr=wsave.getStorage(wsavesave);
2577 : Int ier;
2578 0 : FFTPack::cfft2i(convSize, convSize, wsaveptr, lsav, ier);
2579 : //////////
2580 : #ifdef _OPENMP
2581 0 : omp_set_nested(0);
2582 : #endif
2583 : //////openmp like to share reference param ...making a copy to reduce shared objects
2584 0 : Int cpConvSize=convSize;
2585 0 : Int cpWConvSize=wConvSize;
2586 0 : Double cpWscale=wScale;
2587 0 : Int cpConvSamp=convSampling;
2588 :
2589 : ///////////////
2590 :
2591 0 : convSupport.resize(wConvSize);
2592 0 : convSupport=-1;
2593 0 : Vector<Int> pcsupp;
2594 0 : pcsupp=convSupport;
2595 : Bool delsupstor;
2596 0 : Int* suppstor=pcsupp.getStorage(delsupstor);
2597 : ////////////
2598 :
2599 : //Float max0=1.0;
2600 0 : #pragma omp parallel for default(none) firstprivate(cpWConvSize, cpConvSize, convFuncPtr, s0, s1, wsaveptr, ier, lsav, maxptr, cpWscale,inner, cor, maxConvSize, cpConvSamp, suppstor, C::pi )
2601 :
2602 : for (Int iw=0; iw< cpWConvSize;iw++) {
2603 : // First the w term
2604 : Matrix<Complex> screen(cpConvSize, cpConvSize);
2605 : Matrix<Complex> screen2(cpConvSize, cpConvSize);
2606 : screen=0.0;
2607 : screen2=0.0;
2608 : Bool cpscr;
2609 : Bool cpscr2;
2610 : Complex *scr=screen.getStorage(cpscr);
2611 : Complex *scr2=screen2.getStorage(cpscr2);
2612 : Double twoPiW=2.0*C::pi*Double(iw*iw)/cpWscale;
2613 : for (Int iy=-inner/2;iy<inner/2;iy++) {
2614 : Double m=s1*Double(iy);
2615 : Double msq=m*m;
2616 : //////Int offset= (iy+convSize/2)*convSize;
2617 : ///fftpack likes it flipped
2618 : ooLong offset= (iy>-1 ? iy : ooLong(iy+cpConvSize))*ooLong(cpConvSize);
2619 : for (Int ix=-inner/2;ix<inner/2;ix++) {
2620 : ////// Int ind=offset+ix+convSize/2;
2621 : ///fftpack likes it flipped
2622 : ooLong ind=offset+(ix > -1 ? ooLong(ix) : ooLong(ix+cpConvSize));
2623 : Double l=s0*Double(ix);
2624 : Double rsq=l*l+msq;
2625 : if(rsq<1.0) {
2626 : Double phase=twoPiW*(sqrt(1.0-rsq)-1.0);
2627 : Double cval, sval;
2628 : SINCOS(phase, sval, cval);
2629 :
2630 : Complex comval(cval, sval);
2631 : scr2[ind]=(cor[ooLong(ix+inner/2)+ ooLong((iy+inner/2))*ooLong(inner)])*comval;
2632 : scr[ind]=comval;
2633 :
2634 : }
2635 : }
2636 :
2637 : }
2638 : // Now FFT and get the result back
2639 : /////////Por FFTPack
2640 : Vector<Float>work(2*cpConvSize*cpConvSize);
2641 : Int lenwrk=2*cpConvSize*cpConvSize;
2642 : Bool worksave;
2643 : Float *workptr=work.getStorage(worksave);
2644 : FFTPack::cfft2f(cpConvSize, cpConvSize, cpConvSize, scr, wsaveptr, lsav, workptr, lenwrk, ier);
2645 : FFTPack::cfft2f(cpConvSize, cpConvSize, cpConvSize, scr2, wsaveptr, lsav, workptr, lenwrk, ier);
2646 : screen.putStorage(scr, cpscr);
2647 : screen2.putStorage(scr2, cpscr2);
2648 : ooLong offset=uInt(iw*(cpConvSize/2-1)*(cpConvSize/2-1));
2649 : maxptr[iw]=screen(0,0);
2650 : for (uInt y=0; y< uInt(cpConvSize/2)-1; ++y){
2651 : for (uInt x=0; x< uInt(cpConvSize/2)-1; ++x){
2652 : convFuncPtr[offset+ooLong(y*(cpConvSize/2-1))+ooLong(x)] = screen(x,y);
2653 : }
2654 : }
2655 : /* Bool found=false;
2656 : Int trial=0;
2657 : for (trial=0; trial<cpConvSize/2-2;++trial) {
2658 : // if((abs(convFunc(trial,0,iw))>1e-3)||(abs(convFunc(0,trial,iw))>1e-3) ) {
2659 : if((abs(screen(trial,0))<1e-3)||(abs(screen(0,trial))<1e-3) ) {
2660 : //cout <<"iw " << iw << " x " << abs(convFunc(trial,0,iw)) << " y "
2661 : // <<abs(convFunc(0,trial,iw)) << endl;
2662 : found=true;
2663 : break;
2664 : }
2665 : }
2666 : if(found) {
2667 : suppstor[iw]=Int(0.5+Float(trial)/Float(cpConvSamp))+1;
2668 : if(suppstor[iw]*cpConvSamp*2 >= maxConvSize){
2669 : suppstor[iw]=cpConvSize/2/cpConvSamp-1;
2670 : }
2671 : }
2672 : */
2673 : }
2674 :
2675 :
2676 : #ifdef _OPENMP
2677 0 : omp_set_nested(0);
2678 : #endif
2679 :
2680 : /* convFunc.putStorage(convFuncPtr, convFuncStor);
2681 : maxes.putStorage(maxptr, maxdel);
2682 : Complex maxconv=max(abs(maxes));
2683 : convFunc=convFunc/real(maxconv);
2684 : convFuncPtr=convFunc.getStorage(convFuncStor);
2685 : */
2686 :
2687 0 : #pragma omp parallel for default(none) firstprivate(suppstor, cpConvSize, cpWConvSize, cpConvSamp, convFuncPtr, maxConvSize, maxes)
2688 : for (Int iw=0;iw<cpWConvSize;iw++) {
2689 : Bool found=false;
2690 : Int trial=0;
2691 : ooLong ploffset=ooLong(cpConvSize/2-1)*ooLong(cpConvSize/2-1)*ooLong(iw);
2692 : ////////////////
2693 : // for (trial=cpConvSize/2-2;trial>0;trial--) {
2694 : // // if((abs(convFunc(trial,0,iw))>1e-3)||(abs(convFunc(0,trial,iw))>1e-3) ) {
2695 : // if((abs(convFuncPtr[trial+ploffset])>1e-2)||(abs(convFuncPtr[trial*(cpConvSize/2-1)+ploffset])>1e-2) ) {
2696 : // //cout <<"iw " << iw << " x " << abs(convFunc(trial,0,iw)) << " y "
2697 : // // <<abs(convFunc(0,trial,iw)) << endl;
2698 : // found=true;
2699 : // break;
2700 : // }
2701 : // }
2702 : ///////////////////////
2703 :
2704 : for (trial=0; trial<cpConvSize/2-2;++trial) {
2705 : //for (trial=cpConvSize/2-2;trial>0;trial--) {
2706 : // if((abs(convFunc(trial,0,iw))>1e-3)||(abs(convFunc(0,trial,iw))>1e-3) ) {
2707 : if((abs(convFuncPtr[ooLong(trial)+ploffset])/abs(maxes[0])< 1e-3)||(abs(convFuncPtr[ooLong(trial*(cpConvSize/2-1))+ploffset])/abs(maxes[0]) < 1e-3) ) {
2708 : //if((abs(convFuncPtr[ooLong(trial)+ploffset]) < 1e-3)||(abs(convFuncPtr[ooLong(trial*(cpConvSize/2-1))+ploffset]) < 1e-3) ) {
2709 : //if(abs(convFuncPtr[ooLong(trial)+ploffset])*abs(convFuncPtr[ooLong(trial*(cpConvSize/2-1))+ploffset]) < 1e-3){
2710 : ///diagonal
2711 : //if(abs(convFuncPtr[ooLong((Double(trial)/sqrt(2.0))*(cpConvSize/2-1))+ploffset+ooLong((Double(trial)/sqrt(2.0)))]) < 1e-4){
2712 : //cout <<"iw " << iw << " x " << abs(convFunc(trial,0,iw)) << " y "
2713 : // <<abs(convFunc(0,trial,iw)) << endl;
2714 : found=true;
2715 : break;
2716 : }
2717 : }
2718 :
2719 : if(found) {
2720 : suppstor[iw]=Int(0.5+Float(trial)/Float(cpConvSamp))+1;
2721 : }
2722 : else{
2723 : suppstor[iw]=Int(0.5+Float(cpConvSize/2-2)/Float(cpConvSamp))+1;
2724 : }
2725 : if(suppstor[iw]*cpConvSamp*2 >= maxConvSize){
2726 : suppstor[iw]=cpConvSize/2/cpConvSamp-1;
2727 :
2728 : }
2729 : }
2730 :
2731 :
2732 :
2733 0 : convFunc.putStorage(convFuncPtr, convFuncStor);
2734 : /* maxes.putStorage(maxptr, maxdel);
2735 : Complex maxconv=max(abs(maxes));
2736 : convFunc=convFunc/real(maxconv);
2737 : */
2738 0 : pcsupp.putStorage(suppstor, delsupstor);
2739 0 : convSupport=pcsupp;
2740 : ////////////TESTING
2741 :
2742 : //cerr<< " convsupp " << convSupport << endl;
2743 :
2744 : //convSupport.set(0);
2745 : /////////////////
2746 :
2747 :
2748 : // Normalize such that plane 0 sums to 1 (when jumping in
2749 : // steps of convSampling)
2750 0 : Complex pbSum=0.0;
2751 0 : Int iz=0;
2752 : //for (Int iz=0; iz< convSupport.shape()[0]; ++iz){
2753 0 : pbSum=0.0;
2754 :
2755 0 : for (Int iy=-convSupport(iz);iy<=convSupport(iz);iy++) {
2756 0 : for (Int ix=-convSupport(iz);ix<=convSupport(iz);ix++) {
2757 0 : pbSum+=real(convFunc(abs(ix)*cpConvSamp,abs(iy)*cpConvSamp,iz));
2758 : }
2759 : }
2760 : //convFunc.xyPlane(iz) = convFunc.xyPlane(iz)/Complex(pbSum);
2761 : // }
2762 : //cerr << "pbSum " << pbSum << endl;
2763 0 : for (Int iz2=0; iz2< convSupport.shape()[0]; ++iz2){
2764 0 : convFunc.xyPlane(iz2)= convFunc.xyPlane(iz2)/pbSum;
2765 :
2766 : }
2767 0 : Int newConvSize=2*(max(convSupport)+2)*convSampling;
2768 :
2769 0 : if(newConvSize < convSize){
2770 0 : IPosition blc(3, 0,0,0);
2771 0 : IPosition trc(3, (newConvSize/2-2),
2772 0 : (newConvSize/2-2),
2773 0 : convSupport.shape()(0)-1);
2774 :
2775 0 : Cube<Complex> newConvFunc=convFunc(blc,trc);
2776 0 : convFunc.resize();
2777 0 : convFunc=newConvFunc;
2778 : // convFunctions_p[actualConvIndex_p]->assign(Cube<Complex>(convFunc(blc,trc)));
2779 0 : convSize=newConvSize;
2780 : //cerr << "new convsize " << convSize << endl;
2781 0 : }
2782 :
2783 :
2784 : //////////////////////TESTING
2785 : //convFunc*=Complex(1.0/5.4, 0.0);
2786 :
2787 : ////////////////////////////
2788 : //// if(pbSum>0.0) {
2789 : //// convFunc*=Complex(1.0/pbSum,0.0);
2790 : //// }
2791 : /// else {
2792 : /// cerr << "Convolution function integral is not positive"
2793 : /// << endl;
2794 : /// }
2795 : /////////////TEST
2796 :
2797 : /* Int maxConvSupp=max(convSupport)+1;
2798 : for (Int iw=1;iw<cpWConvSize; ++iw) {
2799 : Float maxplane=max(amplitude(convFunc.xyPlane(iw)));
2800 : for (Int iy=0; iy< maxConvSupp*cpConvSamp; ++iy){
2801 : for (Int ix=0; ix< maxConvSupp*cpConvSamp; ++ix){
2802 : convFunc(ix, iy, iw) /=maxplane;
2803 : }
2804 : }
2805 : }
2806 : */
2807 : /////////////////////////
2808 : /////Write out the SF correction image
2809 : /*String corrim=outMSName_p+String("/WprojCorrection.image");
2810 : Path elpath(corrim);
2811 : cerr << "Saving the correction image " << elpath.absoluteName() << "\nIt should be used to restore images to a flat noise state " << endl;
2812 : if(!Table::isReadable(corrim)){
2813 : ConvolveGridder<Double, Float>elgridder(IPosition(2, nx_p, ny_p),
2814 : uvScale, uvOffset,
2815 : "SF");
2816 : CoordinateSystem newcoord=csys_p;
2817 : StokesCoordinate st(Vector<Int>(1, Stokes::I));
2818 : newcoord.replaceCoordinate(st, 1);
2819 : PagedImage<Float> thisScreen(IPosition(4, nx_p, ny_p, 1, 1), newcoord, corrim);
2820 : thisScreen.set(0.0);
2821 : Vector<Float> correction(nx_p);
2822 : correction=1.0;
2823 :
2824 : Int npixCorr= max(nx_p,ny_p);
2825 : Vector<Float> sincConv(npixCorr);
2826 : for (Int ix=0;ix<npixCorr;ix++) {
2827 : Float x=C::pi*Float(ix-npixCorr/2)/(Float(npixCorr)*Float(convSampling));
2828 : if(ix==npixCorr/2) {
2829 : sincConv(ix)=1.0;
2830 : }
2831 : else {
2832 : sincConv(ix)=sin(x)/x;
2833 : }
2834 : }
2835 : IPosition cursorShape(4, nx_p, 1, 1, 1);
2836 : IPosition axisPath(4, 0, 1, 2, 3);
2837 : LatticeStepper lsx(thisScreen.shape(), cursorShape, axisPath);
2838 : LatticeIterator<Float> lix(thisScreen, lsx);
2839 : for(lix.reset();!lix.atEnd();lix++) {
2840 : Int iy=lix.position()(1);
2841 : elgridder.correctX1D(correction, iy);
2842 :
2843 : for (Int ix=0;ix<nx_p; ++ix) {
2844 : //correction(ix)=sincConv(ix)*sincConv(iy);
2845 : correction(ix)*=sincConv(ix)*sincConv(iy);
2846 : }
2847 : lix.rwVectorCursor()=correction;
2848 : }
2849 : }
2850 :
2851 : */
2852 : // Write out FT of screen as an image
2853 : /* if(1) {
2854 : CoordinateSystem ftCoords(coords);
2855 : Int directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
2856 : AlwaysAssert(directionIndex>=0, AipsError);
2857 : dc=coords.directionCoordinate(directionIndex);
2858 : Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
2859 : Vector<Int> shape(2); shape(0)=convSize;shape(1)=convSize;
2860 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
2861 : ftCoords.replaceCoordinate(*ftdc, directionIndex);
2862 : delete ftdc; ftdc=0;
2863 : ostringstream name;
2864 : name << "FTScreen" ;
2865 : if(Table::canDeleteTable(name)) Table::deleteTable(name);
2866 : PagedImage<Complex> thisScreen(IPosition(4, convFunc.shape()(0), convFunc.shape()(1), 1, convFunc.shape()(2)), ftCoords, name);
2867 : thisScreen.put(convFunc.reform(IPosition(4, convFunc.shape()(0), convFunc.shape()(1), 1, convFunc.shape()(2))));
2868 : if(Table::canDeleteTable("CORR2")) Table::deleteTable("CORR2");
2869 : PagedImage<Complex> thisScreen2(IPosition(4, corr.shape()(0), corr.shape()(1), 1, 1), ftCoords, "CORR2");
2870 : thisScreen2.put(corr.reform(IPosition(4, corr.shape()(0), corr.shape()(1), 1, 1)));
2871 :
2872 : //LatticeExpr<Float> le(real(twoDPB));
2873 : //thisScreen.copyData(le);
2874 : //thisScreen.put(real(screen));
2875 : }
2876 : */
2877 :
2878 :
2879 0 : }
2880 :
2881 0 : Int MSUVBin::sepCommaEmptyToVectorStrings(Vector<String>& lesStrings,
2882 : const String& str){
2883 :
2884 0 : String oneStr=str;
2885 0 : Int nsep=0;
2886 : // decide if its comma seperated or empty space seperated
2887 0 : casacore::String sep;
2888 0 : if((nsep=oneStr.freq(",")) > 0){
2889 0 : sep=",";
2890 : }
2891 : else {
2892 0 : nsep=oneStr.freq(" ");
2893 0 : sep=" ";
2894 : }
2895 0 : if(nsep == 0){
2896 0 : lesStrings.resize(1);
2897 0 : lesStrings=oneStr;
2898 0 : nsep=1;
2899 : }
2900 : else{
2901 0 : String *splitstrings = new String[nsep+1];
2902 0 : nsep=split(oneStr, splitstrings, nsep+1, sep);
2903 0 : lesStrings.resize(nsep);
2904 0 : Int index=0;
2905 0 : for (Int k=0; k < nsep; ++k){
2906 0 : if((String(splitstrings[k]) == String(""))
2907 0 : || (String(splitstrings[k]) == String(" "))){
2908 0 : lesStrings.resize(lesStrings.nelements()-1, true);
2909 : }
2910 : else{
2911 0 : lesStrings[index]=splitstrings[k];
2912 0 : ++index;
2913 : }
2914 : }
2915 0 : delete [] splitstrings;
2916 : }
2917 :
2918 0 : return nsep;
2919 :
2920 0 : }
2921 :
2922 :
2923 :
2924 : } //# NAMESPACE CASA - END
2925 :
|
