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