Line data Source code
1 : //# tAWPFTM.cc: Tests the AWProject class of FTMachines
2 : //# Copyright (C) 2016
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 : #include <casacore/measures/Measures/Stokes.h>
28 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
29 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
30 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
31 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
32 : #include <casacore/coordinates/Coordinates/Projection.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/images/Images/ImageInterface.h>
35 : #include <casacore/images/Images/PagedImage.h>
36 : #include <casacore/images/Images/TempImage.h>
37 : #include <components/ComponentModels/ComponentList.h>
38 : #include <components/ComponentModels/ComponentShape.h>
39 : #include <components/ComponentModels/Flux.h>
40 : #include <casacore/tables/TaQL/ExprNode.h>
41 : #include <casacore/measures/Measures/MeasTable.h>
42 : #include <casacore/ms/MSSel/MSSelection.h>
43 : #include <synthesis/TransformMachines2/SimpleComponentFTMachine.h>
44 : #include <msvis/MSVis/VisImagingWeight.h>
45 : #include <msvis/MSVis/VisibilityIterator2.h>
46 : #include <msvis/MSVis/VisBuffer2.h>
47 : #include <casacore/casa/namespace.h>
48 : #include <casacore/casa/OS/Directory.h>
49 : #include <casacore/casa/Utilities/Regex.h>
50 : #include <synthesis/TransformMachines/FTMachine.h>
51 : #include <synthesis/TransformMachines2/FTMachine.h>
52 : #include <synthesis/TransformMachines2/AWProjectWBFTNew.h>
53 : #include <synthesis/TransformMachines2/AWProjectWBFTHPG.h>
54 : //#include <synthesis/TransformMachines2/AWConvFuncEPJones.h>
55 : #include <synthesis/TransformMachines2/EVLAAperture.h>
56 : #include <synthesis/TransformMachines2/AWVisResampler.h>
57 : #include <synthesis/TransformMachines2/AWVisResamplerHPG.h>
58 : #include <synthesis/TransformMachines2/PointingOffsets.h>
59 : #include <synthesis/TransformMachines2/VBStore.h>
60 : #include <casacore/images/Images/ImageUtilities.h>
61 : #include <casacore/images/Images/ImageOpener.h>
62 : #include <casacore/images/Images/ImageSummary.h>
63 : #include <casacore/casa/System/ProgressMeter.h>
64 : #include <casacore/casa/OS/Timer.h>
65 : #include <msvis/MSVis/ViFrequencySelection.h>
66 : #include <casacore/casa/OS/EnvVar.h>
67 : #include <cstdlib>
68 : //#include <synthesis/TransformMachines2/HPGModelImage.h>
69 :
70 : using namespace casa;
71 : using namespace casacore;
72 : using namespace casa::refim;
73 : using namespace casacore;
74 : using namespace std;
75 : //using namespace casa::test;
76 : /*
77 : void createAWPFTMachine(const String ftmName,
78 : const String modelImageName,
79 : CountedPtr<refim::FTMachine>& theFT, CountedPtr<refim::FTMachine>& theIFT,
80 : const String& telescopeName,
81 : MPosition& observatoryLocation,
82 : const String cfCache= "testCF.cf",
83 : const Int cfBufferSize=1024,
84 : const Int cfOversampling=20,
85 : const Bool wbAWP= true,
86 : //------------------------------
87 : const Int wprojPlane=1,
88 : //const Float=1.0,
89 : const Bool useDoublePrec=true,
90 : //------------------------------
91 : const Bool aTermOn= true,
92 : const Bool psTermOn= false,
93 : const Bool mTermOn= false,
94 : const Bool doPointing= true,
95 : const Bool doPBCorr= true,
96 : const Bool conjBeams= true,
97 : Float pbLimit_l=1e-3,
98 : vector<float> posigdev = {300.0,300.0},
99 : const Float computePAStep=360.0,
100 : const Float rotatePAStep=360.0,
101 : const Int cache=1000000000,
102 : const Int tile=16,
103 : const String imageNamePrefix=String(""))
104 :
105 : {
106 : LogIO os( LogOrigin("tAWPFTM","createAWPFTMachine",WHERE));
107 :
108 : if (wprojPlane<=1)
109 : {
110 : os << LogIO::NORMAL
111 : << "You are using wprojplanes=1. Doing co-planar imaging (no w-projection needed)"
112 : << LogIO::POST;
113 : os << LogIO::NORMAL << "Performing WBA-Projection" << LogIO::POST; // Loglevel PROGRESS
114 : }
115 : else
116 : os << LogIO::NORMAL << "Performing WBAW-Projection" << LogIO::POST; // Loglevel PROGRESS
117 : // if((wprojPlane>1)&&(wprojPlane<64))
118 : // {
119 : // os << LogIO::WARN
120 : // << "No. of w-planes set too low for W projection - recommend at least 128"
121 : // << LogIO::POST;
122 : // }
123 :
124 : // MSObservationColumns msoc(mss4vi_p[0].observation());
125 : // String telescopeName=msoc.telescopeName()(0);
126 : CountedPtr<refim::ConvolutionFunction> awConvFunc = AWProjectFT::makeCFObject(telescopeName,
127 : aTermOn,
128 : psTermOn, true, mTermOn, wbAWP,
129 : false);
130 : CountedPtr<refim::PointingOffsets> po = new refim::PointingOffsets();
131 : po->setOverSampling(awConvFunc->getOversampling());
132 : po->setDoPointing(doPointing);
133 : // awConvFunc->setPointingOffsets(po);
134 : //
135 : // Construct the appropriate re-sampler.
136 : //
137 : CountedPtr<refim::VisibilityResamplerBase> visResampler;
138 : if (ftmName == "awphpg")
139 : visResampler = new refim::AWVisResamplerHPG();
140 : else
141 : visResampler = new refim::AWVisResampler();
142 :
143 : visResampler->setModelImage(modelImageName);
144 : //
145 : // Construct and initialize the CF cache object.
146 : //
147 : CountedPtr<refim::CFCache> cfCacheObj;
148 :
149 : //
150 : // Finally construct the FTMachine with the CFCache, ConvFunc and
151 : // Re-sampler objects.
152 : //
153 : // Float pbLimit_l=1e-3;
154 : // vector<float> posigdev = {300.0,300.0};
155 : // Float pbLimit_l=1e-3;
156 : if(ftmName=="awphpg"){
157 : theFT=new refim::AWProjectWBFTHPG(wprojPlane, cache/2,
158 : cfCacheObj, awConvFunc,
159 : visResampler,
160 : doPointing, posigdev, doPBCorr,
161 : tile, computePAStep, pbLimit_l, true,conjBeams,
162 : useDoublePrec);
163 : theFT->setPBReady(true);
164 : theIFT=theFT;
165 : }
166 : else
167 : {
168 : theFT = new refim::AWProjectWBFTNew(wprojPlane, cache/2,
169 : cfCacheObj, awConvFunc,
170 : visResampler,
171 : true */ /*doPointing, posigdev, doPBCorr,
172 : tile, computePAStep, pbLimit_l, true,conjBeams,
173 : useDoublePrec);
174 : }
175 : cfCacheObj = new refim::CFCache();
176 : cfCacheObj->setCacheDir(cfCache.data()); /////TESTOO LAZY FILL ?
177 : cfCacheObj->setLazyFill(False);
178 :
179 : cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
180 : cfCacheObj->initCache2();
181 :
182 : theFT->setCFCache(cfCacheObj);
183 : */
184 : /*
185 : Quantity rotateOTF(rotatePAStep,"deg");
186 : static_cast<refim::AWProjectWBFTNew &>(*theFT).setObservatoryLocation(observatoryLocation);
187 : static_cast<refim::AWProjectWBFTNew &>(*theFT).setPAIncrement(Quantity(computePAStep,"deg"),rotateOTF);
188 :
189 : theIFT = new refim::AWProjectWBFTNew(static_cast<refim::AWProjectWBFTNew &>(*theFT));
190 : //// Send in Freq info.
191 : // os << "Sending frequency selection information " << mssFreqSel_p << " to AWP FTM." << LogIO::POST;
192 : // theFT->setSpwFreqSelection( mssFreqSel_p );
193 : // theIFT->setSpwFreqSelection( mssFreqSel_p );
194 :
195 : }
196 :
197 : //
198 :
199 : bool mdFromString(casacore::MDirection &theDir, const casacore::String &in)
200 : {
201 : bool rstat(false);
202 : String tmpA, tmpB, tmpC;
203 : std::istringstream iss(in);
204 : iss >> tmpA >> tmpB >> tmpC;
205 : casacore::Quantity tmpQA;
206 : casacore::Quantity tmpQB;
207 : casacore::Quantity::read(tmpQA, tmpA);
208 : casacore::Quantity::read(tmpQB, tmpB);
209 : if(tmpC.length() > 0){
210 : MDirection::Types theRF;
211 : MDirection::getType(theRF, tmpC);
212 : theDir = MDirection (tmpQA, tmpQB, theRF);
213 : rstat = true;
214 : } else {
215 : theDir = MDirection (tmpQA, tmpQB);
216 : rstat = true;
217 : }
218 : return rstat;
219 : }
220 :
221 : std::map<casacore::Int, std::map<casacore::Int, casacore::Vector<casacore::Int> > > makeTheChanSelMap(MSSelection& mss)
222 : {
223 : Matrix<Int> mssChanSel = mss.getChanList();
224 : IPosition mssChanSelShape=mssChanSel.shape();
225 :
226 : std::map<casacore::Int, std::map<casacore::Int, casacore::Vector<casacore::Int> > > chanselKurfuffel;
227 : std::map<casacore::Int, casacore::Vector<casacore::Int> > channelsPerSpw;
228 : for(int i=0;i<mssChanSelShape[0];i++)
229 : {
230 : Vector<Int> channels(2,-1);
231 : int spw=mssChanSel(i,0);
232 : channels[0]=mssChanSel(i,2)-mssChanSel(i,1)+1;//mssChanSel(i,1);
233 : channels[1]=mssChanSel(i,1);//mssChanSel(i,2);
234 : channelsPerSpw[spw]=channels;
235 : // cerr << "SPW:CHAN " << spw << " " << channels << endl;
236 : }
237 : chanselKurfuffel[0]=channelsPerSpw;
238 : return chanselKurfuffel;
239 : }
240 : PagedImage<Complex> makeEmptySkyImage(VisibilityIterator2& vi2,
241 : const MeasurementSet& selectedMS,
242 : MSSelection& msSelection,
243 : const String& imageName, const String& startModelImageName,
244 : const Vector<Int>& imSize, const float& cellSize,
245 : const String& phaseCenter, const String& stokes,
246 : const String& refFreq,
247 : const String& mode)
248 : {
249 : Vector<Quantity> qCellSize(2);
250 : for (int i=0;i<2;i++)
251 : {
252 : String qUnit("arcsec");
253 : qCellSize[i]=Quantity(cellSize, qUnit);
254 : }
255 :
256 : int imStokes=1;
257 : if (stokes=="I") imStokes=1;
258 : else if (stokes=="IV") imStokes=2;
259 : else if (stokes=="IQUV") imStokes=4;
260 :
261 : int imNChan=1;
262 : if (mode=="mfs") imNChan=1;
263 : if(mode=="cube") imNChan=3;
264 : // else if (mode=="pseudo") {}
265 : // else if (mode=="spectral") {imnchan=datanchan[0];imstart=datastart[0];imstep=datastep[0];}
266 :
267 : MDirection mphaseCenter;
268 : mdFromString(mphaseCenter, phaseCenter);
269 :
270 : SynthesisParamsGrid gridParams;
271 : SynthesisParamsImage imageParams;
272 : imageParams.imageName=imageName;
273 : imageParams.startModel=startModelImageName;
274 : imageParams.imsize=imSize;
275 : imageParams.cellsize=qCellSize;
276 : imageParams.phaseCenter = mphaseCenter;
277 : imageParams.stokes=stokes;
278 : imageParams.mode=mode;
279 : imageParams.frame=String("LSRK");
280 : imageParams.veltype=String("radio");
281 : if(mode=="cube"){
282 : imageParams.start="1.2GHz";
283 : imageParams.step="300MHz";
284 : imageParams.freqStart=Quantity(1.2, "GHz");
285 : imageParams.freqStep=Quantity(300, "MHz");
286 : }
287 :
288 : //
289 : // There are two items related to ref. freq. "reffreq" (a String) and "refFreq" (a Quantity).
290 : // Set them both.
291 : //
292 : if (refFreq != "")
293 : {
294 : imageParams.reffreq=refFreq;
295 :
296 : std::istringstream iss(refFreq);
297 : Double ff;
298 : iss >> ff;
299 : // cerr << "Ref. freq. = " << ff << endl;
300 : imageParams.refFreq=Quantity(ff/1e9,"GHz");
301 : }
302 :
303 : casacore::Block<const casacore::MeasurementSet *> msList(1); msList[0]=&selectedMS;
304 : CoordinateSystem csys = imageParams.buildCoordinateSystem(vi2,makeTheChanSelMap(msSelection),msList);
305 : IPosition imshape(4,imSize(0),imSize(1),imStokes,imNChan);
306 : return PagedImage<Complex>(imshape, csys, imageParams.imageName);
307 : }
308 : //
309 : //-------------------------------------------------------------------------
310 : //
311 : std::tuple<Vector<Int>, Vector<Int> > loadMS(const String& msname,
312 : const String& spwStr,
313 : const String& fieldStr,
314 : const String& uvDistStr,
315 : MeasurementSet& thems,
316 : MeasurementSet& selectedMS,
317 : MSSelection& msSelection)
318 : {
319 : //MeasurementSet thems;
320 : if (Table::isReadable(msname))
321 : thems=MeasurementSet(msname, Table::Update);
322 : else
323 : throw(AipsError(msname+" does exist or not readable"));
324 : //
325 : //-------------------------------------------------------------------
326 : //
327 : msSelection.setSpwExpr(spwStr);
328 : //msSelection.setAntennaExpr(antStr);
329 : msSelection.setFieldExpr(fieldStr);
330 : msSelection.setUvDistExpr(uvDistStr);
331 : selectedMS = MeasurementSet(thems);
332 : Vector<int> spwid, fieldid;
333 : TableExprNode exprNode=msSelection.toTableExprNode(&thems);
334 : if (!exprNode.isNull())
335 : {
336 : selectedMS = MS(thems(exprNode));
337 : spwid=msSelection.getSpwList();
338 : fieldid=msSelection.getFieldList();
339 : }
340 : return std::tuple<Vector<Int>, Vector<Int> >{spwid, fieldid};
341 : }
342 : //
343 : //-------------------------------------------------------------------------
344 : //
345 : */
346 1 : Int main(int argc, char **argv)
347 : {
348 : /*
349 : String casadata = EnvironmentVariable::get("CASADATA");
350 : if (casadata.empty()) {
351 : throw(AipsError(
352 : "CASADATA env variable not defined and no file given in command line"));
353 : }
354 : //initialize beamCalc with repo path
355 : BeamCalc::Instance(String(casadata+"/../distro/"));
356 : String source_input_ms =
357 : casadata + "/measurementset/evla/refim_mawproject.ms";
358 : //
359 : // -------------------------------------- Just the UI -------------------------------------------------------------------
360 : //
361 : system("rm -rf ./refim_mawproject.ms");
362 :
363 : system(String("cp -r "+source_input_ms+" .").c_str());
364 : string MSNBuf="./refim_mawproject.ms";
365 : string ftmName="awphpg";
366 : string cfCache="test.cf";
367 : string fieldStr="";
368 : string spwStr="*";
369 : string uvDistStr="";
370 : string imageName="gooloo";
371 : string modelImageName="";
372 : string cmplxGridName="gooloo.cmplx";
373 : string phaseCenter="19h59m28.5 40d44m01.5 J2000";
374 : string stokes="I";
375 : string refFreqStr="3.0e9";
376 : string weighting="natural";
377 :
378 : float refFreq=3e09, freqBW=3e9;
379 : float cellSize=10;
380 : int NX=512;
381 : Int NY=512;
382 : Int cfBufferSize=512;
383 : Int cfOversampling=20;
384 : Int nW=1;
385 : int nInt=200;
386 : int nchan=1;
387 : Bool WBAwp=true;
388 : Bool restartUI=false;
389 : Bool doPointing=false;
390 : Bool normalize=false;
391 : Bool doPBCorr= false;
392 : Bool conjBeams= false;
393 : Float pbLimit=1e-3;
394 : vector<float> posigdev = {300.0,300.0};
395 :
396 : //
397 : // -------------------------------------- End of UI -------------------------------------------------------------------
398 : //
399 :
400 : // cerr << "###Info: Pointing sigma dev = " << posigdev[0] << "," << //posigdev[1] << endl;
401 : // restartUI=false;
402 : NY=NX;
403 :
404 : std::stringstream strNChan,strNInt;
405 :
406 : strNChan << nchan;
407 : strNInt << nInt;
408 : Timer timer;
409 : double griddingTime=0;
410 : try
411 : {
412 : if (phaseCenter=="") throw(AipsError("The phasecenter parameter needs to be set"));
413 : if (refFreqStr=="") throw(AipsError("The reffreq parameter needs to be set"));
414 : if (NX<=0) throw(AipsError("The imsize parameter needs to be set to a positive finite value"));
415 : if (cellSize<=0) throw(AipsError("The cellsize parameter needs to be set to a positive finite value"));
416 :
417 : {
418 : Directory dirObj(modelImageName);
419 : if ((modelImageName!="") && ((!dirObj.exists()) || (!dirObj.isReadable())))
420 : {
421 : std::string msg;
422 : msg="Model image \""+modelImageName + "\" not found.";
423 : throw(AipsError(msg));
424 : }
425 : }
426 :
427 : //-------------------------------------------------------------------
428 : // Load the selected MS. The original ms (thems), the selected
429 : // MS and the MSSelection objects are modified. The selected
430 : // list of SPW and FIELD IDs are returned as a std::tuple.
431 : //
432 : hpg::initialize();
433 : MSSelection msSelection;
434 : MS theMS, selectedMS;
435 : {
436 : // String msname=MSNBuf;
437 : Vector<int> spwid, fieldid;
438 : LogIO log_l(LogOrigin("tAWPHPG","main"));
439 :
440 : log_l << "Opening the MS (\"" << MSNBuf << "\"), applying data selection, "
441 : << "setting up data iterators...all that boring stuff."
442 : << LogIO::POST;
443 : auto lists = loadMS(MSNBuf,
444 : spwStr,
445 : fieldStr,
446 : uvDistStr,
447 : theMS,
448 : selectedMS,
449 : msSelection);
450 : spwid=std::get<0>(lists);
451 : fieldid=std::get<1>(lists);
452 : }
453 : //
454 : //-------------------------------------------------------------------
455 : //
456 :
457 : //-------------------------------------------------------------------
458 : // Set up the data iterator
459 : //
460 : vi::VisibilityIterator2 vi2(selectedMS,vi::SortColumns(),true,0,10.0);
461 : {
462 : Matrix<Double> freqSelection= msSelection.getChanFreqList(NULL,true);
463 : // vi2.setInterval(50.0);
464 : VisImagingWeight viw(weighting);
465 : vi2.useImagingWeight(viw);
466 : FrequencySelectionUsingFrame fsel(MFrequency::Types::LSRK);
467 : for(unsigned i=0;i<freqSelection.shape()(0); i++)
468 : fsel.add(freqSelection(i,0), freqSelection(i,1), freqSelection(i,2));
469 : vi2.setFrequencySelection (fsel);
470 : }
471 : vi::VisBuffer2 *vb=vi2.getVisBuffer();
472 :
473 : //
474 :
475 :
476 :
477 : //-------------------------------------------------------------------
478 : // Make the empty grid with the sky image coordinates
479 : //
480 : Vector<Int> imSize(2,NX);
481 : String mode="mfs", startModelImageName="";
482 :
483 : // String gridName="cgrid.tmp";
484 : // if (cmplxGridName != "") gridName=cmplxGridName;
485 :
486 : casacore::MDirection pc;
487 : Int pcField = casa::refim::getPhaseCenter(selectedMS,pc);
488 : std::ostringstream oss;
489 : pc.print(oss);
490 : // cerr << "PC = " << oss << endl;
491 : ComponentList cl;
492 : SkyComponent otherPoint(ComponentType::POINT);
493 : otherPoint.flux() = Flux<Double>(6.66e-2, 0.0, 0.0, 0.00000);
494 : otherPoint.shape().setRefDirection(pc);
495 : cl.add(otherPoint);
496 : PagedImage<Complex> cgrid=makeEmptySkyImage(vi2, selectedMS, msSelection,
497 : cmplxGridName, startModelImageName,
498 : imSize, cellSize, phaseCenter,
499 : stokes, refFreqStr, mode);
500 : PagedImage<Float> skyImage(cgrid.shape(),cgrid.coordinates(), imageName);
501 : PagedImage<Float> weightImage(cgrid.shape(),cgrid.coordinates(), imageName+".weight");
502 : if (cmplxGridName=="")
503 : cgrid.table().markForDelete();
504 :
505 : StokesImageUtil::From(cgrid, skyImage);
506 : if(vb->polarizationFrame()==MSIter::Linear) StokesImageUtil::changeCStokesRep(cgrid,StokesImageUtil::LINEAR);
507 : else StokesImageUtil::changeCStokesRep(cgrid, StokesImageUtil::CIRCULAR);
508 :
509 : //
510 : //-------------------------------------------------------------------
511 : //
512 :
513 : //-------------------------------------------------------------------
514 : // Create the AWP FTMachine. The AWProjectionFT is construed
515 : // with the re-sampler depending on the ftmName (AWVisResampler
516 : // for ftmName='awproject' or AWVisResamplerHPG for
517 : // ftmName='awphpg')
518 : //
519 : CountedPtr<refim::FTMachine> ftm,iftm;
520 : //FTMachine * ftm=nullptr;
521 : MPosition loc;
522 : MeasTable::Observatory(loc, MSColumns(selectedMS).observation().telescopeName()(0));
523 : Bool useDoublePrec=true, aTermOn=true, psTermOn=false, mTermOn=false;
524 :
525 : createAWPFTMachine(ftmName, modelImageName, ftm, iftm,
526 : String("EVLA"),
527 : loc, cfCache, cfBufferSize, cfOversampling,WBAwp,nW,
528 : useDoublePrec,
529 : aTermOn,
530 : psTermOn,
531 : mTermOn,
532 : doPointing,
533 : doPBCorr,
534 : conjBeams,
535 : pbLimit,
536 : posigdev
537 : );
538 : {
539 : Matrix<Double> mssFreqSel;
540 : mssFreqSel = msSelection.getChanFreqList(NULL,true);
541 : // Send in Freq info.
542 : // cerr << "Sending frequency selection information " << mssFreqSel << " to AWP FTM." << endl;
543 : ftm->setSpwFreqSelection( mssFreqSel );
544 : iftm->setSpwFreqSelection( mssFreqSel );
545 : }
546 : //
547 : //-------------------------------------------------------------------
548 : //
549 :
550 : //-------------------------------------------------------------------
551 : // Iterative over the MS, trigger the gridding.
552 : //
553 : {
554 : cgrid.set(Complex(0.0));
555 : Matrix<Float> weight;
556 : vi2.originChunks();
557 : vi2.origin();
558 : ftm->initializeToSky(cgrid, weight, *vb);
559 : //drygridding
560 : /* {
561 : cerr << "@@@@@DRY GRIDDING" << endl;
562 : ftm->setDryRun(true);
563 : for (vi2.originChunks();vi2.moreChunks(); vi2.nextChunk())
564 : {
565 : for (vi2.origin(); vi2.more(); vi2.next()){
566 : ftm->put(*vb, -1, True, casa::refim::FTMachine::PSF);
567 : }
568 :
569 : }
570 : Directory dir(cfCache);
571 : Vector<String> cfList=dir.find(Regex(Regex::fromPattern("CFS*")));
572 : cerr << "CFLIST " << cfList << endl;
573 : Vector<String> wtCFList;
574 : wtCFList.resize(cfList.nelements());
575 : for (Int i=0; i<(Int)wtCFList.nelements(); i++) wtCFList[i]="WT"+cfList[i];
576 :
577 : auto cfc = ftm->getCFCache();
578 : cfc->initCacheFromList2(cfCache, cfList, wtCFList,
579 : 360.0, 360.0,False);
580 :
581 : Vector<Double> uvScale, uvOffset;
582 : Matrix<Double> vbFreqSelection;
583 : CountedPtr<refim::CFStore2> cfs2 = CountedPtr<refim::CFStore2>(&cfc->memCache2_p[0],false);
584 : CountedPtr<refim::CFStore2> cfwts2 = CountedPtr<refim::CFStore2>(&cfc->memCacheWt2_p[0],false);
585 : casa::refim::AWConvFunc::makeConvFunction2(cfCache, uvScale, uvOffset, vbFreqSelection,
586 : *cfs2, *cfwts2, True, True, False);
587 : ftm->setDryRun(false);
588 : }
589 : */
590 : // cerr << "image.shape: " << cgrid.shape() << endl;
591 : /*
592 : refim::SimpleComponentFTMachine cft;
593 :
594 : timer.mark();
595 : int n=0;
596 : //Timer tvi;
597 : // Vector<refim::VBStore> vbsList;
598 : // vbsList.resize(10);
599 : unsigned long vol=0;
600 : ProgressMeter pm(1.0, selectedMS.nrow(),
601 : "Gridding", "","","",true);
602 :
603 : casa::refim::FTMachine::Type dataCol=casa::refim::FTMachine::CORRECTED;
604 : if(selectedMS.tableDesc().isColumn("CORRECTED_DATA") ) dataCol = casa::refim::FTMachine::CORRECTED;
605 : else
606 : {
607 : LogIO log_l(LogOrigin("tAWPHPG","main"));
608 : log_l << "CORRECTED_DATA column not found. Using the DATA column instead." << LogIO::WARN << LogIO::POST;
609 : dataCol = casa::refim::FTMachine::OBSERVED;
610 : }
611 :
612 : for (vi2.originChunks();vi2.moreChunks(); vi2.nextChunk())
613 : {
614 : for (vi2.origin(); vi2.more(); vi2.next())
615 : {
616 : //tvi.mark();
617 :
618 : //casacore::Cube<casacore::Complex> data(vb->visCube());
619 : //vol+=data.shape().product()*sizeof(casacore::Complex);
620 :
621 : //cft.get(*vb, cl);
622 : //vb->setVisCube(vb->visCubeModel());
623 :
624 : // Check if CORRECTED_DATA column exits. Hope is that
625 : // the VisBuffer does reasonable caching (i.e., the
626 : // hope is that the code below does not read the data
627 : // twice).
628 : if (dataCol==casa::refim::FTMachine::CORRECTED)
629 : vb->setVisCube(vb->visCubeCorrected());
630 : else
631 : vb->setVisCube(vb->visCube());
632 : ftm->setFTMType(dataCol);
633 : ftm->put(*vb, -1, False, dataCol);
634 :
635 : vol+=vb->nRows();
636 : pm.update(Double(vol));
637 : n++;
638 : // cerr << "Iter " << n << " " << tvi.real() << endl;
639 : // if (n > 0) break;
640 : }
641 : // if (n > 0) break;
642 : }
643 : griddingTime += timer.real();
644 : cerr << "Gridding time: " << griddingTime << ". No. of rows processed: " << vol << ". Data rate: " << vol/griddingTime << " rows/s" << endl;
645 : ftm->finalizeToSky();
646 :
647 : // Get the weights. Do nothing with them for now.
648 : //Bool normalize=false;
649 : ftm->getImage(weight, normalize);
650 : // Convert the skyImage (retrived in ftm->finalizeToSky()) and
651 : // convert it from Feed basis to Stokes basis.
652 : StokesImageUtil::To(skyImage, cgrid);
653 :
654 : cerr << "val at center " << cgrid.getAt(IPosition(4, NX/2, NY/2, 0, 0)) << " sumweight " << weight << endl;
655 :
656 : //Lets do the weight image
657 : cgrid.set(Complex(0.0));
658 : weight.resize();
659 : vi2.originChunks();
660 : vi2.origin();
661 : vol=0;
662 : ftm->initializeToSky(cgrid, weight, *vb);
663 : for (vi2.originChunks();vi2.moreChunks(); vi2.nextChunk())
664 : {
665 : for (vi2.origin(); vi2.more(); vi2.next())
666 : {
667 :
668 :
669 : ftm->setFTMType(refim::FTMachine::WEIGHT);
670 : ftm->put(*vb, -1, False);
671 :
672 : vol+=vb->nRows();
673 : pm.update(Double(vol));
674 : //n++;
675 : // cerr << "Iter " << n << " " << tvi.real() << endl;
676 : // if (n > 0) break;
677 : }
678 : // if (n > 0) break;
679 :
680 : }
681 :
682 : ftm->finalizeToSky();
683 :
684 : // Get the weights. Do nothing with them for now.
685 : //Bool normalize=false;
686 : ftm->getWeightImage(weightImage, weight);
687 : StokesImageUtil::To(skyImage, cgrid);
688 :
689 : cerr << "val at center " << weightImage.getAt(IPosition(4, NX/2, NY/2, 0, 0)) << " sumweight " << weight << endl;
690 :
691 : }
692 : ///////
693 : //detach the ms for cleaning up
694 : selectedMS=MeasurementSet();
695 : }
696 : catch (AipsError& x)
697 : {
698 : cerr << "AipsError exception: " << x.getMesg() << endl;
699 : // RestartUI(RENTER);
700 : return(1);
701 : }
702 : catch(...)
703 : {
704 : throw(AipsError("Some kind of Exception caught"));
705 : //x << "###AipsError: " << "\t\"" << x.what() << "\"" << endl;
706 : // RestartUI(RENTER);
707 : return(1);
708 : }
709 : cerr <<"OK" << endl;
710 : hpg::finalize();
711 : exit(0);
712 : */
713 1 : }
714 :
|
