Commits

Mark Kettenis authored c64451fed80 Merge
Merge branch 'master' into CAS-13270

casa5/code/msvis/MSVis/VisImagingWeight.cc

Modified
90 90 String mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId());
91 91 multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(mapid, 0));
92 92 gwt_p[0]=new TempLattice<Float>(IPosition(2, nx,ny), 0);
93 93 gwt_p[0]->set(0.0);
94 94 a_gwt_p.resize(nx_p, ny_p);
95 95 a_gwt_p.set(0.0);
96 96
97 97 Int fields=0;
98 98 for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
99 99 for (vi.origin();vi.more();vi++) {
100 -
100 +
101 101 if(vb->newFieldId()){
102 102 mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId());
103 103 if(multiField){
104 104 if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) ){
105 105 fields+=1;
106 106 gwt_p.resize(fields+1);
107 107 gwt_p[fields]=new TempLattice<Float>(IPosition(2, nx_p,ny_p),0);
108 108 gwt_p[fields]->set(0.0);
109 109 }
110 110 }
115 115 }
116 116 }
117 117
118 118 Float u, v;
119 119 Vector<Double> sumwt(fields+1,0.0);
120 120 f2_p.resize(fields+1);
121 121 d2_p.resize(fields+1);
122 122 Int fid=0;
123 123 for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
124 124 for (vi.origin();vi.more();vi++) {
125 - if(vb->newFieldId()){
125 + if(vb->newFieldId()){
126 126 mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId());
127 127
128 128
129 - if(multiField){
130 - if(activeFieldIndex_p >=0)
131 - gwt_p[activeFieldIndex_p]->put(a_gwt_p);
132 - }
133 - activeFieldIndex_p=multiFieldMap_p.find(mapid) !=multiFieldMap_p.end( ) ? multiFieldMap_p[mapid] : -1;
134 - if(multiField && activeFieldIndex_p >=0)
135 - gwt_p[activeFieldIndex_p]->get(a_gwt_p);
136 -
137 - }
129 + if(multiField){
130 + if(activeFieldIndex_p >=0)
131 + gwt_p[activeFieldIndex_p]->put(a_gwt_p);
132 + }
133 + activeFieldIndex_p=multiFieldMap_p.find(mapid) !=multiFieldMap_p.end( ) ? multiFieldMap_p[mapid] : -1;
134 + if(multiField && activeFieldIndex_p >=0)
135 + gwt_p[activeFieldIndex_p]->get(a_gwt_p);
136 +
137 + }
138 138 auto mapvp = multiFieldMap_p.find(mapid);
139 139 fid= mapvp != multiFieldMap_p.end( ) ? mapvp->second : -1;
140 140 Int nRow=vb->nRow();
141 141 Int nChan=vb->nChannel();
142 142
143 - // Extract weights correctly (use WEIGHT_SPECTRUM, if available)
144 - Matrix<Float> wtm;
145 - ////////////////
146 - doWtSp=(vb->weightSpectrum().nelements()>0);
147 - //////////////
148 - if (doWtSp)
149 - // WS available,
150 - unPolChanWeight(wtm,vb->weightSpectrum()); // collapse corr axis (parallel-hand average)
151 - else
152 - // WS UNavailable
153 - wtm.reference(vb->weight().reform(IPosition(2,1,nRow))); // use vb.weight() (corr-collapsed, w/ 1 channel)
154 -
155 - // Use this to mod the channel indexing below
156 - Int nChanWt=wtm.shape()(0);
143 + // Extract weights correctly (use WEIGHT_SPECTRUM, if available)
144 + Matrix<Float> wtm;
145 + ////////////////
146 + doWtSp=(vb->weightSpectrum().nelements()>0);
147 + //////////////
148 + if (doWtSp)
149 + // WS available,
150 + unPolChanWeight(wtm,vb->weightSpectrum()); // collapse corr axis (parallel-hand average)
151 + else
152 + // WS UNavailable
153 + wtm.reference(vb->weight().reform(IPosition(2,1,nRow))); // use vb.weight() (corr-collapsed, w/ 1 channel)
154 +
155 + // Use this to mod the channel indexing below
156 + Int nChanWt=wtm.shape()(0);
157 157
158 158 for (Int row=0; row<nRow; row++) {
159 159 for (Int chn=0; chn<nChan; chn++) {
160 160 if(!vb->flag()(chn,row)) {
161 - Float currwt=wtm(chn%nChanWt,row); // the weight for this chan,row
161 + Float currwt=wtm(chn%nChanWt,row); // the weight for this chan,row
162 162 Float f=vb->frequency()(chn)/C::c;
163 163 u=vb->uvw()(row)(0)*f;
164 164 v=vb->uvw()(row)(1)*f;
165 165 Int ucell=Int(uscale_p*u+uorigin_p);
166 166 Int vcell=Int(vscale_p*v+vorigin_p);
167 -
167 +
168 168 if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) {
169 169 for (Int iv=-vBox;iv<=vBox;iv++) {
170 170 for (Int iu=-uBox;iu<=uBox;iu++) {
171 - a_gwt_p(ucell+iu,vcell+iv)+=currwt;
171 + a_gwt_p(ucell+iu,vcell+iv)+=currwt;
172 172 sumwt[fid]+=currwt;
173 173 }
174 174 }
175 175 }
176 176 ucell=Int(-uscale_p*u+uorigin_p);
177 177 vcell=Int(-vscale_p*v+vorigin_p);
178 178 if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) {
179 179 for (Int iv=-vBox;iv<=vBox;iv++) {
180 180 for (Int iu=-uBox;iu<=uBox;iu++) {
181 181 a_gwt_p(ucell+iu,vcell+iv)+=currwt;
191 191 // make sure last active plane is saved
192 192 gwt_p[activeFieldIndex_p]->put(a_gwt_p);
193 193 // We use the approximation that all statistical weights are equal to
194 194 // calculate the average summed weights (over visibilities, not bins!)
195 195 // This is simply to try an ensure that the normalization of the robustness
196 196 // parameter is similar to that of the ungridded case, but it doesn't have
197 197 // to be exact, since any given case will require some experimentation.
198 198
199 199 //Float f2, d2;
200 200 for(fid=0; fid < Int(gwt_p.nelements()); ++fid){
201 - gwt_p[fid]->get(a_gwt_p);
202 - activeFieldIndex_p=fid;
203 - if (rmode=="norm") {
201 + gwt_p[fid]->get(a_gwt_p);
202 + activeFieldIndex_p=fid;
203 + if (rmode=="norm") {
204 204 os << "Normal robustness, robust = " << robust << LogIO::POST;
205 205 Double sumlocwt = 0.;
206 206 for(Int vgrid=0;vgrid<ny;vgrid++) {
207 207 for(Int ugrid=0;ugrid<nx;ugrid++) {
208 208 if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
209 209 }
210 210 }
211 211 f2_p[fid] = square(5.0*pow(10.0,Double(-robust))) / (sumlocwt / sumwt[fid]);
212 212 d2_p[fid] = 1.0;
213 213
230 230 VisImagingWeight::VisImagingWeight(ROVisibilityIterator& vi, Block<CountedPtr<TempLattice<Float> > >& grids, const String& rmode, const Quantity& noise,
231 231 const Double robust, const Quantity& cellx, const Quantity& celly,
232 232 const Bool multiField) : doFilter_p(false), robust_p(robust), rmode_p(rmode), noise_p(noise) {
233 233
234 234 LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE));
235 235
236 236 VisBufferAutoPtr vb (vi);
237 237 wgtType_p="uniform";
238 238 gwt_p.resize(grids.nelements(), true, false);
239 239 for (uInt k =0; k < gwt_p.nelements(); ++k)
240 - gwt_p[k]=grids[k];
240 + gwt_p[k]=grids[k];
241 241 nx_p=gwt_p[0]->shape()[0];
242 242 ny_p=gwt_p[0]->shape()[1];
243 243 uscale_p=(nx_p*cellx.get("rad").getValue());
244 244 vscale_p=(ny_p*celly.get("rad").getValue());
245 245 uorigin_p=nx_p/2;
246 246 vorigin_p=ny_p/2;
247 247
248 248 multiFieldMap_p.clear();
249 249 Int fields=0;
250 250 String mapid;
251 251 for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
252 - for (vi.origin();vi.more();vi++) {
253 - if(vb->newFieldId()){
254 - mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId());
255 - if(multiField){
256 - if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) ){
257 - fields+=1;
252 + for (vi.origin();vi.more();vi++) {
253 + if(vb->newFieldId()){
254 + mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId());
255 + if(multiField){
256 + if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) ){
257 + fields+=1;
258 258
259 - }
260 - }
261 - if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) )
262 - multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(mapid, fields));
263 - }
264 - }
259 + }
260 + }
261 + if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) )
262 + multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(mapid, fields));
263 + }
264 + }
265 265 }
266 266 f2_p.resize(gwt_p.nelements());
267 267 d2_p.resize(gwt_p.nelements());
268 268 for(uInt fid=0; fid < gwt_p.nelements(); ++fid){
269 - activeFieldIndex_p=fid;
270 - gwt_p[fid]->get(a_gwt_p);
271 - if (rmode_p=="norm") {
272 - os << "Normal robustness, robust = " << robust_p << LogIO::POST;
273 - Double sumlocwt = 0.;
274 - for(Int vgrid=0;vgrid<ny_p;vgrid++) {
275 - for(Int ugrid=0;ugrid<nx_p;ugrid++) {
276 - if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
277 - }
278 - }
279 - Double sumwt_fid=sum(a_gwt_p);
280 - f2_p[fid] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumwt_fid);
281 - d2_p[fid] = 1.0;
282 -
283 - }
269 + activeFieldIndex_p=fid;
270 + gwt_p[fid]->get(a_gwt_p);
271 + if (rmode_p=="norm") {
272 + os << "Normal robustness, robust = " << robust_p << LogIO::POST;
273 + Double sumlocwt = 0.;
274 + for(Int vgrid=0;vgrid<ny_p;vgrid++) {
275 + for(Int ugrid=0;ugrid<nx_p;ugrid++) {
276 + if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
277 + }
278 + }
279 + Double sumwt_fid=sum(a_gwt_p);
280 + f2_p[fid] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumwt_fid);
281 + d2_p[fid] = 1.0;
282 +
283 + }
284 284 else if (rmode=="abs") {
285 285 os << "Absolute robustness, robust = " << robust_p << ", noise = "
286 286 << noise_p.get("Jy").getValue() << "Jy" << LogIO::POST;
287 287 f2_p[fid] = square(robust_p);
288 288 d2_p[fid] = 2.0 * square(noise_p.get("Jy").getValue());
289 289
290 290 }
291 291 else {
292 292 f2_p[fid] = 1.0;
293 293 d2_p[fid] = 0.0;
294 294 }
295 295 }
296 296 }
297 297
298 298 VisImagingWeight::VisImagingWeight(vi::VisibilityIterator2& visIter,
299 - const String& rmode, const Quantity& noise,
299 + const String& rmode, const Quantity& noise,
300 300 const Double robust, const Int nx, const Int ny,
301 301 const Quantity& cellx, const Quantity& celly,
302 302 const Int uBox, const Int vBox, const Bool multiField) : doFilter_p(false), robust_p(robust), rmode_p(rmode), noise_p(noise), activeFieldIndex_p(-1) {
303 303
304 304 LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE));
305 305
306 306
307 307
308 308 vi::VisBuffer2* vb=(visIter.getVisBuffer());
309 309 wgtType_p="uniform";
353 353 }
354 354 }
355 355
356 356 Float u, v;
357 357 Vector<Double> sumwt(fields+1,0.0);
358 358 f2_p.resize(fields+1);
359 359 d2_p.resize(fields+1);
360 360 Int fid=0;
361 361 for (visIter.originChunks();visIter.moreChunks();visIter.nextChunk()) {
362 362 for (visIter.origin();visIter.more();visIter.next()) {
363 - if(vb->isNewFieldId()){
363 + if(vb->isNewFieldId()){
364 364 mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]);
365 365
366 - if(multiField){
367 - if(activeFieldIndex_p >=0)
368 - gwt_p[activeFieldIndex_p]->put(a_gwt_p);
369 - }
370 - activeFieldIndex_p= multiFieldMap_p.find(mapid) != multiFieldMap_p.end( ) ? multiFieldMap_p[mapid] : -1;
371 - if(multiField && activeFieldIndex_p >=0)
372 - gwt_p[activeFieldIndex_p]->get(a_gwt_p);
373 - }
366 + if(multiField){
367 + if(activeFieldIndex_p >=0)
368 + gwt_p[activeFieldIndex_p]->put(a_gwt_p);
369 + }
370 + activeFieldIndex_p= multiFieldMap_p.find(mapid) != multiFieldMap_p.end( ) ? multiFieldMap_p[mapid] : -1;
371 + if(multiField && activeFieldIndex_p >=0)
372 + gwt_p[activeFieldIndex_p]->get(a_gwt_p);
373 + }
374 374
375 375 auto mapvp = multiFieldMap_p.find(mapid);
376 376 fid= mapvp != multiFieldMap_p.end( ) ? mapvp->second : -1;
377 377 Int nRow=vb->nRows();
378 378 Int nChan=vb->nChannels();
379 379
380 - // Extract weights correctly (use WEIGHT_SPECTRUM, if available)
381 - Matrix<Float> wtm;
382 - Cube<Float> wtc;
383 - if (doWtSp)
384 - // WS available [nCorr,nChan,nRow]
385 - wtc.reference(vb->weightSpectrum());
386 - else {
387 - Int nCorr=vb->nCorrelations();
388 - // WS UNavailable weight()[nCorr,nRow] --> [nCorr,nChan,nRow]
389 - wtc.reference(vb->weight().reform(IPosition(3,nCorr,1,nRow))); // unchan'd weight as single-chan
390 - }
391 - unPolChanWeight(wtm,wtc); // Collapse on corr axis
392 -
393 - // Used for mod of chn index below
394 - Int nChanWt=wtm.shape()(0);
395 -
396 - //Oww !!! temporary implementation of old vb.flag just to see if things work
397 - Matrix<Bool> flag;
398 - cube2Matrix(vb->flagCube(), flag);
380 + // Extract weights correctly (use WEIGHT_SPECTRUM, if available)
381 + Matrix<Float> wtm;
382 + Cube<Float> wtc;
383 + if (doWtSp)
384 + // WS available [nCorr,nChan,nRow]
385 + wtc.reference(vb->weightSpectrum());
386 + else {
387 + Int nCorr=vb->nCorrelations();
388 + // WS UNavailable weight()[nCorr,nRow] --> [nCorr,nChan,nRow]
389 + wtc.reference(vb->weight().reform(IPosition(3,nCorr,1,nRow))); // unchan'd weight as single-chan
390 + }
391 + unPolChanWeight(wtm,wtc); // Collapse on corr axis
392 +
393 + // Used for mod of chn index below
394 + Int nChanWt=wtm.shape()(0);
395 +
396 + //Oww !!! temporary implementation of old vb.flag just to see if things work
397 + Matrix<Bool> flag;
398 + cube2Matrix(vb->flagCube(), flag);
399 399 for (Int row=0; row<nRow; row++) {
400 400 for (Int chn=0; chn<nChan; chn++) {
401 -
402 - Float currwt=wtm(chn%nChanWt,row); // the weight for this chan,row
403 -
401 +
402 + Float currwt=wtm(chn%nChanWt,row); // the weight for this chan,row
403 +
404 404 if(!flag(chn,row)) {
405 405 Float f=vb->getFrequency(row, chn)/C::c;
406 406 u=vb->uvw()(0,row)*f;
407 407 v=vb->uvw()(1,row)*f;
408 408 Int ucell=Int(uscale_p*u+uorigin_p);
409 409 Int vcell=Int(vscale_p*v+vorigin_p);
410 410 if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) {
411 411 for (Int iv=-vBox;iv<=vBox;iv++) {
412 412 for (Int iu=-uBox;iu<=uBox;iu++) {
413 - a_gwt_p(ucell+iu,vcell+iv)+=currwt;
413 + a_gwt_p(ucell+iu,vcell+iv)+=currwt;
414 414 sumwt[fid]+=currwt;
415 415 }
416 416 }
417 417 }
418 418 ucell=Int(-uscale_p*u+uorigin_p);
419 419 vcell=Int(-vscale_p*v+vorigin_p);
420 420 if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) {
421 421 for (Int iv=-vBox;iv<=vBox;iv++) {
422 422 for (Int iu=-uBox;iu<=uBox;iu++) {
423 423 a_gwt_p(ucell+iu,vcell+iv)+=currwt;
434 434 // make sure last active plane is saved
435 435 gwt_p[activeFieldIndex_p]->put(a_gwt_p);
436 436 // We use the approximation that all statistical weights are equal to
437 437 // calculate the average summed weights (over visibilities, not bins!)
438 438 // This is simply to try an ensure that the normalization of the robustness
439 439 // parameter is similar to that of the ungridded case, but it doesn't have
440 440 // to be exact, since any given case will require some experimentation.
441 441
442 442 //Float f2, d2;
443 443 for(fid=0; fid < Int(gwt_p.nelements()); ++fid){
444 - gwt_p[fid]->get(a_gwt_p);
445 - activeFieldIndex_p=fid;
446 - if (rmode=="norm") {
444 + gwt_p[fid]->get(a_gwt_p);
445 + activeFieldIndex_p=fid;
446 + if (rmode=="norm") {
447 447 os << "Normal robustness, robust = " << robust << LogIO::POST;
448 448 Double sumlocwt = 0.;
449 449 for(Int vgrid=0;vgrid<ny;vgrid++) {
450 450 for(Int ugrid=0;ugrid<nx;ugrid++) {
451 451 if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
452 452 }
453 453 }
454 454 f2_p[fid] = square(5.0*pow(10.0,Double(-robust))) / (sumlocwt / sumwt[fid]);
455 455 d2_p[fid] = 1.0;
456 456
483 483 dc.setWorldAxisUnits(Vector<String>(2, "rad"));
484 484 Double cellx=dc.increment()(0);
485 485 Double celly=dc.increment()(1);
486 486 uscale_p=nx_p*cellx;
487 487 vscale_p=ny_p*celly;
488 488 uorigin_p=nx_p/2;
489 489 vorigin_p=ny_p/2;
490 490 //Now to recover from image density and other parameters
491 491 Int nplanes=1;
492 492 if(im.shape().nelements()==5)
493 - nplanes=im.shape()[4];
493 + nplanes=im.shape()[4];
494 494 gwt_p.resize(nplanes, True, False);
495 495 if(im.shape().nelements()==5){
496 - IPosition blc(Vector<Int>(5,0));
497 - for (Int fid=0;fid<nplanes;fid++)
498 - {
499 - gwt_p[fid]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
500 - Array<Float> lala;
501 - blc[4]=fid;
502 - im.getSlice(lala, blc, IPosition(5, nx_p, ny_p,1,1,1), True);
503 - gwt_p[fid]->put( lala.reform(IPosition(2, nx_p, ny_p)));
504 -
505 - }
506 - }
507 - else{
508 - Array<Float> lala;
509 - im.get(lala, True);
510 - gwt_p[0]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
511 - gwt_p[0]->put( lala.reform(IPosition(2, nx_p, ny_p)));
512 -
513 - }
496 + IPosition blc(Vector<Int>(5,0));
497 + for (Int fid=0;fid<nplanes;fid++)
498 + {
499 + gwt_p[fid]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
500 + Array<Float> lala;
501 + blc[4]=fid;
502 + im.getSlice(lala, blc, IPosition(5, nx_p, ny_p,1,1,1), True);
503 + gwt_p[fid]->put( lala.reform(IPosition(2, nx_p, ny_p)));
504 +
505 + }
506 + }
507 + else{
508 + Array<Float> lala;
509 + im.get(lala, True);
510 + gwt_p[0]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
511 + gwt_p[0]->put( lala.reform(IPosition(2, nx_p, ny_p)));
512 +
513 + }
514 514 const TableRecord& rec=im.miscInfo();
515 515 if(rec.isDefined("d2")){
516 - d2_p.resize();
517 - rec.get("d2", d2_p);
518 - f2_p.resize();
519 - rec.get("f2", f2_p);
520 - multiFieldMap_p.clear();
516 + d2_p.resize();
517 + rec.get("d2", d2_p);
518 + f2_p.resize();
519 + rec.get("f2", f2_p);
520 + multiFieldMap_p.clear();
521 521 Int mapsize;
522 522 rec.get("multimapsize", mapsize);
523 - for(Int k=0; k < mapsize; ++k){
524 - String key;
525 - Int val;
526 - rec.get("key"+String::toString(k), key);
527 - rec.get("val"+String::toString(k), val);
523 + for(Int k=0; k < mapsize; ++k){
524 + String key;
525 + Int val;
526 + rec.get("key"+String::toString(k), key);
527 + rec.get("val"+String::toString(k), val);
528 528 //cerr << "key and id " << key << " "<< val << endl;
529 - multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(key, val));
530 - }
531 -
529 + multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(key, val));
530 + }
531 +
532 532
533 533 }
534 534 activeFieldIndex_p=0;
535 535 gwt_p[0]->get(a_gwt_p);
536 536 if(rec.isDefined("dofilter")){
537 537 rec.get("dofilter", doFilter_p);
538 538 rec.get("rbmaj", rbmaj_p);
539 539 rec.get("rbmin", rbmin_p);
540 540 rec.get("cospa", cospa_p);
541 541 rec.get("sinpa", sinpa_p);
542 542 }
543 543
544 544
545 545 }
546 546
547 547
548 548 VisImagingWeight::~VisImagingWeight(){
549 549 /*for (uInt fid=0; fid < gwt_p.nelements(); ++fid){
550 550 gwt_p[fid].resize();
551 - }*/
551 + }*/
552 552
553 553 }
554 554
555 555 Vector<Int> VisImagingWeight::shapeOfdensityGrid(){
556 556 Vector<Int> retval(3, 0);
557 557 retval(2)=gwt_p.nelements();
558 558 if(retval(2) > 0){
559 - retval[0]=gwt_p[0]->shape()(0);
560 - retval[1]=gwt_p[0]->shape()(1);
559 + retval[0]=gwt_p[0]->shape()(0);
560 + retval[1]=gwt_p[0]->shape()(1);
561 561
562 562 }
563 563
564 564 return retval;
565 565 }
566 566 void VisImagingWeight::toImageInterface(casacore::ImageInterface<casacore::Float>& im){
567 567
568 568 if( wgtType_p != "uniform")
569 - throw(AipsError("cannot save weight density for non-Briggs' weighting schemes"));
569 + throw(AipsError("cannot save weight density for non-Briggs' weighting schemes"));
570 570
571 - IPosition where= IPosition(im.shape().nelements(),0);
571 + IPosition where= IPosition(im.shape().nelements(),0);
572 572 Int lastAx=where.nelements()-1;
573 573 for (uInt fid=0;fid<gwt_p.nelements(); ++fid){
574 - activeFieldIndex_p=fid;
575 - gwt_p[fid]->get(a_gwt_p);
576 - where[lastAx]=fid;
577 - im.putSlice(a_gwt_p, where);
574 + activeFieldIndex_p=fid;
575 + gwt_p[fid]->get(a_gwt_p);
576 + where[lastAx]=fid;
577 + im.putSlice(a_gwt_p, where);
578 578
579 579 }
580 580 Record rec(im.miscInfo());
581 581 rec.define("d2", d2_p);
582 582 rec.define("f2", f2_p);
583 583 rec.define("numfield", Int(multiFieldMap_p.size()));
584 584 uInt keycount=0;
585 585 for (auto iter = multiFieldMap_p.begin( ); iter != multiFieldMap_p.end( ); ++iter, ++keycount){
586 - rec.define("key"+String::toString(keycount), iter->first);
587 - rec.define("val"+String::toString(keycount), iter->second);
586 + rec.define("key"+String::toString(keycount), iter->first);
587 + rec.define("val"+String::toString(keycount), iter->second);
588 588 }
589 589 rec.define("multimapsize",keycount);
590 590 if(doFilter_p){
591 591 rec.define("dofilter", doFilter_p);
592 592 rec.define("rbmaj", rbmaj_p);
593 593 rec.define("rbmin", rbmin_p);
594 594 rec.define("cospa", cospa_p);
595 595 rec.define("sinpa", sinpa_p);
596 596 }
597 597
598 598
599 599 im.setMiscInfo(rec);
600 600
601 601 }
602 602 void VisImagingWeight::setFilter(const String& type, const Quantity& bmaj,
603 - const Quantity& bmin, const Quantity& bpa)
603 + const Quantity& bmin, const Quantity& bpa)
604 604 {
605 605
606 606 LogIO os(LogOrigin("VisImagingWeight", "setFilter()", WHERE));
607 607
608 608 if (type=="gaussian") {
609 609
610 610 Bool lambdafilt=false;
611 611
612 612 if( bmaj.getUnit().contains("lambda"))
613 - lambdafilt=true;
613 + lambdafilt=true;
614 614 if(lambdafilt){
615 - os << "Filtering for Gaussian of shape: "
616 - << bmaj.get("klambda").getValue() << " by "
617 - << bmin.get("klambda").getValue() << " (klambda) at p.a. "
618 - << bpa.get("deg").getValue() << " (degrees)" << LogIO::POST;
619 - rbmaj_p=log(2.0)/square(bmaj.get("lambda").getValue());
620 - rbmin_p=log(2.0)/square(bmin.get("lambda").getValue());
615 + os << "Filtering for Gaussian of shape: "
616 + << bmaj.get("klambda").getValue() << " by "
617 + << bmin.get("klambda").getValue() << " (klambda) at p.a. "
618 + << bpa.get("deg").getValue() << " (degrees)" << LogIO::POST;
619 + rbmaj_p=log(2.0)/square(bmaj.get("lambda").getValue());
620 + rbmin_p=log(2.0)/square(bmin.get("lambda").getValue());
621 621 }
622 622 else{
623 - os << "Filtering for Gaussian of shape: "
624 - << bmaj.get("arcsec").getValue() << " by "
625 - << bmin.get("arcsec").getValue() << " (arcsec) at p.a. "
626 - << bpa.get("deg").getValue() << " (degrees)" << LogIO::POST;
627 -
628 - // Convert to values that we can use
629 - Double fact = 4.0*log(2.0);
630 - rbmaj_p = fact*square(bmaj.get("rad").getValue());
631 - rbmin_p = fact*square(bmin.get("rad").getValue());
623 + os << "Filtering for Gaussian of shape: "
624 + << bmaj.get("arcsec").getValue() << " by "
625 + << bmin.get("arcsec").getValue() << " (arcsec) at p.a. "
626 + << bpa.get("deg").getValue() << " (degrees)" << LogIO::POST;
627 +
628 + // Convert to values that we can use
629 + Double fact = 4.0*log(2.0);
630 + rbmaj_p = fact*square(bmaj.get("rad").getValue());
631 + rbmin_p = fact*square(bmin.get("rad").getValue());
632 632 }
633 633 Double rbpa = MVAngle(bpa).get("rad").getValue();
634 634 cospa_p = sin(rbpa);
635 635 sinpa_p = cos(rbpa);
636 636 doFilter_p=true;
637 637
638 638 }
639 639 else {
640 - os << "Unknown filtering " << type << LogIO::EXCEPTION;
640 + os << "Unknown filtering " << type << LogIO::EXCEPTION;
641 641 }
642 642
643 643
644 644
645 645
646 646 }
647 647
648 648
649 649 Bool VisImagingWeight::doFilter() const{
650 650
651 651 return doFilter_p;
652 652 }
653 653
654 654
655 - void VisImagingWeight::filter(Matrix<Float>& imWeight, const Matrix<Bool>& flag,
656 - const Matrix<Double>& uvw,
657 - const Vector<Double>& frequency, const Matrix<Float>& weight) const{
655 + void VisImagingWeight::filter(Matrix<Float>& imWeight, const Matrix<Bool>& flag,
656 + const Matrix<Double>& uvw,
657 + const Vector<Double>& frequency, const Matrix<Float>& weight) const{
658 658
659 659 // Expecting weight[nchan,nrow], where nchan=1 or nChan (data)
660 660 // If nchan=1 (i.e., WEIGHT_SPECTRUM unavailable), then the
661 661 // following will be handle the channel degeneracy correctly, using %.
662 662
663 663
664 664 Int nRow=imWeight.shape()(1);
665 665 Int nChan=imWeight.shape()(0);
666 666 Int nChanWt=weight.shape()(0);
667 667 for (Int row=0; row<nRow; row++) {
668 668 for (Int chn=0; chn<nChan; chn++) {
669 - Double invLambdaC=frequency(chn)/C::c;
670 - Double u = uvw(0,row);
671 - Double v = uvw(1,row);
672 - if(!flag(chn,row) && (weight(chn%nChanWt,row)>0.0) ) {
673 - Double ru = invLambdaC*( cospa_p * u + sinpa_p * v);
674 - Double rv = invLambdaC*(- sinpa_p * u + cospa_p * v);
675 - Double filter = exp(-rbmaj_p*square(ru) - rbmin_p*square(rv));
676 - imWeight(chn,row)*=filter;
677 - }
678 - else {
679 - imWeight(chn,row)=0.0;
680 - }
669 + Double invLambdaC=frequency(chn)/C::c;
670 + Double u = uvw(0,row);
671 + Double v = uvw(1,row);
672 + if(!flag(chn,row) && (weight(chn%nChanWt,row)>0.0) ) {
673 + Double ru = invLambdaC*( cospa_p * u + sinpa_p * v);
674 + Double rv = invLambdaC*(- sinpa_p * u + cospa_p * v);
675 + Double filter = exp(-rbmaj_p*square(ru) - rbmin_p*square(rv));
676 + imWeight(chn,row)*=filter;
677 + }
678 + else {
679 + imWeight(chn,row)=0.0;
680 + }
681 681 }
682 682 }
683 683
684 684
685 685 }
686 686
687 687
688 688 void VisImagingWeight::weightUniform(Matrix<Float>& imWeight, const Matrix<Bool>& flag, const Matrix<Double>& uvw,
689 689 const Vector<Double>& frequency,
690 690 const Matrix<Float>& weight, const Int msId, const Int fieldId) const{
691 691
692 692
693 693 //cout << " WEIG " << nx_p << " " << ny_p << " " << gwt_p[0].shape() << endl;
694 - //cout << "f2 " << f2_p[0] << " d2 " << d2_p[0] << " uscale " << uscale_p << " vscal " << vscale_p << " origs " << uorigin_p << " " << vorigin_p << endl;
694 + //cout << "f2 " << f2_p[0] << " d2 " << d2_p[0] << " uscale " << uscale_p << " vscal " << vscale_p << " origs " << uorigin_p << " " << vorigin_p << endl;
695 695 String mapid=String::toString(msId)+String("_")+String::toString(fieldId);
696 696 //cout << "min max gwt " << min(gwt_p[0]) << " " << max(gwt_p[0]) << " mapid " << mapid <<endl;
697 697
698 698 if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) )
699 - throw(AipsError("Imaging weight calculation is requested for a data that was not selected"));
699 + throw(AipsError("Imaging weight calculation is requested for a data that was not selected"));
700 700
701 701 auto mapvp = multiFieldMap_p.find(mapid);
702 702 Int fid = mapvp != multiFieldMap_p.end( ) ? mapvp->second : -1;
703 703 //Int ndrop=0;
704 704 if(activeFieldIndex_p != fid){
705 - a_gwt_p=gwt_p[fid]->get();
706 - activeFieldIndex_p=fid;
705 + a_gwt_p=gwt_p[fid]->get();
706 + activeFieldIndex_p=fid;
707 707 }
708 708 Double sumwt=0.0;
709 709 Int nRow=imWeight.shape()(1);
710 710 Int nChannel=imWeight.shape()(0);
711 711 Int nChanWt=weight.shape()(0);
712 712 Float u, v;
713 713 for (Int row=0; row<nRow; row++) {
714 - for (Int chn=0; chn<nChannel; chn++) {
715 - if (!flag(chn,row)) {
716 - Float f=frequency(chn)/C::c;
717 - u=uvw(0, row)*f;
718 - v=uvw(1, row)*f;
719 - Int ucell=Int(uscale_p*u+uorigin_p);
720 - Int vcell=Int(vscale_p*v+vorigin_p);
721 - imWeight(chn,row)=weight(chn%nChanWt,row);
722 - if((ucell>0)&&(ucell<nx_p)&&(vcell>0)&&(vcell<ny_p) &&a_gwt_p(ucell,vcell)>0.0) {
723 - imWeight(chn,row)/=a_gwt_p(ucell,vcell)*f2_p[fid]+d2_p[fid];
724 - sumwt+=imWeight(chn,row);
725 -
726 - }
727 - else {
728 - imWeight(chn,row)=0.0;
729 - //ndrop++;
730 - }
731 - }
732 - else{
733 - imWeight(chn,row)=0.0;
734 - }
735 - }
714 + for (Int chn=0; chn<nChannel; chn++) {
715 + if (!flag(chn,row)) {
716 + Float f=frequency(chn)/C::c;
717 + u=uvw(0, row)*f;
718 + v=uvw(1, row)*f;
719 + Int ucell=Int(uscale_p*u+uorigin_p);
720 + Int vcell=Int(vscale_p*v+vorigin_p);
721 + imWeight(chn,row)=weight(chn%nChanWt,row);
722 + if((ucell>0)&&(ucell<nx_p)&&(vcell>0)&&(vcell<ny_p) &&a_gwt_p(ucell,vcell)>0.0) {
723 + imWeight(chn,row)/=a_gwt_p(ucell,vcell)*f2_p[fid]+d2_p[fid];
724 + sumwt+=imWeight(chn,row);
725 +
726 + }
727 + else {
728 + imWeight(chn,row)=0.0;
729 + //ndrop++;
730 + }
731 + }
732 + else{
733 + imWeight(chn,row)=0.0;
734 + }
735 + }
736 736 }
737 737
738 738 }
739 739
740 740 /* unused version?
741 741 void VisImagingWeight::weightNatural(Matrix<Float>& imagingWeight, const Matrix<Bool>& flag,
742 742 const Matrix<Float>& weight) const{
743 743
744 744
745 745
746 746 Int nRow=imagingWeight.shape()(1);
747 747 Vector<Float> wgtRow(nRow);
748 -
748 +
749 749 for (Int row=0; row<nRow; row++) {
750 - wgtRow(row)=max(weight.column(row));
750 + wgtRow(row)=max(weight.column(row));
751 751 }
752 - weightNatural(imagingWeight, flag, wgtRow);
752 + weightNatural(imagingWeight, flag, wgtRow);
753 753
754 754 }
755 755 */
756 756 void VisImagingWeight::weightNatural(Matrix<Float>& imagingWeight, const Matrix<Bool>& flag,
757 757 const Matrix<Float>& weight) const{
758 758
759 759 Double sumwt=0.0;
760 - //cerr << "shape of weight " << weight.shape() << " max " << max (weight.column(0) ) << " max " << min(weight.column(0)) << " " << weight.column(0).shape() << endl;
760 + //cerr << "shape of weight " << weight.shape() << " max " << max (weight.column(0) ) << " max " << min(weight.column(0)) << " " << weight.column(0).shape() << endl;
761 761 Int nRow=imagingWeight.shape()(1);
762 762 Int nChan=imagingWeight.shape()(0);
763 763 Int nChanWt=weight.shape()(0);
764 764 for (Int row=0; row<nRow; row++) {
765 765 for (Int chn=0; chn<nChan; chn++) {
766 766 if( !flag(chn,row) ) {
767 - imagingWeight(chn,row)=weight(chn%nChanWt,row);
767 + imagingWeight(chn,row)=weight(chn%nChanWt,row);
768 768 sumwt+=imagingWeight(chn,row);
769 769 }
770 770 else {
771 771 imagingWeight(chn,row)=0.0;
772 772 }
773 773 }
774 774 }
775 775
776 776
777 777 }
778 778
779 779
780 780 void VisImagingWeight::weightRadial(Matrix<Float>& imagingWeight,
781 781 const Matrix<Bool>& flag,
782 782 const Matrix<Double>& uvw,
783 783 const Vector<Double>& frequency,
784 784 const Matrix<Float>& weight) const{
785 785
786 786 Double sumwt=0.0;
787 787 Int nRow=imagingWeight.shape()(1);
788 788 Int nChan=imagingWeight.shape()(0);
789 - Int nChanWt=weight.shape()(0);
789 + Int nChanWt=weight.shape()(0);
790 790
791 791 for (Int row=0; row<nRow; row++) {
792 792 for (Int chn=0; chn< nChan; chn++) {
793 793 Float f=frequency(chn)/C::c;
794 794 if( !flag(chn,row) ) {
795 795 imagingWeight(chn,row)=
796 - f*sqrt(square(uvw(0, row))+square(uvw(1, row)))
797 - *weight(chn%nChanWt,row);
796 + f*sqrt(square(uvw(0, row))+square(uvw(1, row)))
797 + *weight(chn%nChanWt,row);
798 798 sumwt+=imagingWeight(chn,row);
799 799 }
800 800 else {
801 801 imagingWeight(chn,row)=0.0;
802 802 }
803 803 }
804 804 }
805 805
806 806
807 807 }
808 808
809 809 VisImagingWeight& VisImagingWeight::operator=(const VisImagingWeight& other){
810 810 if(this != &other){
811 - // gwt_p=other.gwt_p;
811 + // gwt_p=other.gwt_p;
812 812
813 - gwt_p.resize(other.gwt_p.nelements(), true, false);
814 - for (uInt k=0; k < gwt_p.nelements(); ++k){
815 - //gwt_p[k].resize();
816 - gwt_p[k]=other.gwt_p[k];
817 - }
813 + gwt_p.resize(other.gwt_p.nelements(), true, false);
814 + for (uInt k=0; k < gwt_p.nelements(); ++k){
815 + //gwt_p[k].resize();
816 + gwt_p[k]=other.gwt_p[k];
817 + }
818 818
819 819
820 820 wgtType_p=other.wgtType_p;
821 821 uscale_p=other.uscale_p;
822 822 vscale_p=other.vscale_p;
823 - f2_p.resize();
824 - d2_p.resize();
823 + f2_p.resize();
824 + d2_p.resize();
825 825 f2_p=other.f2_p;
826 826 d2_p=other.d2_p;
827 827 uorigin_p=other.uorigin_p;
828 828 vorigin_p=other.vorigin_p;
829 829 nx_p=other.nx_p;
830 830 ny_p=other.ny_p;
831 - doFilter_p=other.doFilter_p;
832 - cospa_p=other.cospa_p;
833 - sinpa_p=other.sinpa_p;
834 - rbmaj_p=other.rbmaj_p;
835 - rbmin_p=other.rbmin_p;
836 - robust_p=other.robust_p;
837 - rmode_p=other.rmode_p;
838 - multiFieldMap_p=other.multiFieldMap_p;
831 + doFilter_p=other.doFilter_p;
832 + cospa_p=other.cospa_p;
833 + sinpa_p=other.sinpa_p;
834 + rbmaj_p=other.rbmaj_p;
835 + rbmin_p=other.rbmin_p;
836 + robust_p=other.robust_p;
837 + rmode_p=other.rmode_p;
838 + multiFieldMap_p=other.multiFieldMap_p;
839 839 }
840 840 return *this;
841 841 }
842 842
843 843 String VisImagingWeight::getType() const{
844 844
845 845 return wgtType_p;
846 846
847 847 }
848 848 Bool VisImagingWeight::getWeightDensity (Block<Matrix<Float> >& density){
854 854 for (uInt k=0; k < gwt_p.nelements(); ++k){
855 855 density[k].resize(gwt_p[k]->shape());
856 856 density[k]=gwt_p[k]->get();
857 857 }
858 858 return true;
859 859 }
860 860 void VisImagingWeight::setWeightDensity(const Block<Matrix<Float> >& density){
861 861 if(wgtType_p=="uniform"){
862 862 gwt_p.resize(density.nelements(), true, false);
863 863 for (uInt k=0; k < gwt_p.nelements(); ++k){
864 - //gwt_p[k].resize();
865 - gwt_p[k]=new TempLattice<Float>(IPosition(2, density[k].shape()[0], density[k].shape()[1]),0);
866 - gwt_p[k]->put(density[k]);
864 + //gwt_p[k].resize();
865 + gwt_p[k]=new TempLattice<Float>(IPosition(2, density[k].shape()[0], density[k].shape()[1]),0);
866 + gwt_p[k]->put(density[k]);
867 867 }
868 868 //break any old reference
869 869 f2_p.resize();
870 870 d2_p.resize();
871 871 f2_p.resize(gwt_p.nelements());
872 872 d2_p.resize(gwt_p.nelements());
873 873 f2_p.set(1.0);
874 874 d2_p.set(0.0);
875 875 //Float f2, d2;
876 876 for(uInt fid=0; fid < gwt_p.nelements(); ++fid){
877 - activeFieldIndex_p=fid;
878 - a_gwt_p=gwt_p[fid]->get();
879 - if (rmode_p=="norm") {
880 - Double sumlocwt = 0.;
881 - for(Int vgrid=0;vgrid<a_gwt_p.shape()(1);vgrid++) {
882 - for(Int ugrid=0;ugrid<a_gwt_p.shape()(0);ugrid++) {
883 - if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
884 - }
885 - }
886 - Double sumwt_fid=sum(a_gwt_p);
887 - f2_p[fid] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumwt_fid);
888 - d2_p[fid] = 1.0;
889 - }
890 - else if (rmode_p=="abs") {
891 - f2_p[fid] = square(robust_p);
892 - d2_p[fid] = 2.0 * square(noise_p.get("Jy").getValue());
893 -
894 - }
895 - else {
896 - f2_p[fid] = 1.0;
897 - d2_p[fid] = 0.0;
898 - }
877 + activeFieldIndex_p=fid;
878 + a_gwt_p=gwt_p[fid]->get();
879 + if (rmode_p=="norm") {
880 + Double sumlocwt = 0.;
881 + for(Int vgrid=0;vgrid<a_gwt_p.shape()(1);vgrid++) {
882 + for(Int ugrid=0;ugrid<a_gwt_p.shape()(0);ugrid++) {
883 + if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
884 + }
885 + }
886 + Double sumwt_fid=sum(a_gwt_p);
887 + f2_p[fid] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumwt_fid);
888 + d2_p[fid] = 1.0;
889 + }
890 + else if (rmode_p=="abs") {
891 + f2_p[fid] = square(robust_p);
892 + d2_p[fid] = 2.0 * square(noise_p.get("Jy").getValue());
893 +
894 + }
895 + else {
896 + f2_p[fid] = 1.0;
897 + d2_p[fid] = 0.0;
898 + }
899 899 }
900 900
901 901 }
902 902
903 903
904 904 }
905 905 void VisImagingWeight::cube2Matrix(const Cube<Bool>& fcube, Matrix<Bool>& fMat)
906 906 {
907 - fMat.resize(fcube.shape()[1], fcube.shape()[2]);
908 - Bool deleteIt1;
909 - Bool deleteIt2;
910 - const Bool * pcube = fcube.getStorage (deleteIt1);
911 - Bool * pflags = fMat.getStorage (deleteIt2);
912 - for (uInt row = 0; row < fcube.shape()[2]; row++) {
913 - for (Int chn = 0; chn < fcube.shape()[1]; chn++) {
914 - *pflags = *pcube++;
915 - for (Int pol = 1; pol < fcube.shape()[0]; pol++, pcube++) {
916 - *pflags = *pcube ? *pcube : *pflags;
917 - }
918 - pflags++;
919 - }
920 - }
921 - fcube.freeStorage (pcube, deleteIt1);
922 - fMat.putStorage (pflags, deleteIt2);
907 + fMat.resize(fcube.shape()[1], fcube.shape()[2]);
908 + Bool deleteIt1;
909 + Bool deleteIt2;
910 + const Bool * pcube = fcube.getStorage (deleteIt1);
911 + Bool * pflags = fMat.getStorage (deleteIt2);
912 + for (uInt row = 0; row < fcube.shape()[2]; row++) {
913 + for (Int chn = 0; chn < fcube.shape()[1]; chn++) {
914 + *pflags = *pcube++;
915 + for (Int pol = 1; pol < fcube.shape()[0]; pol++, pcube++) {
916 + *pflags = *pcube ? *pcube : *pflags;
917 + }
918 + pflags++;
919 + }
920 + }
921 + fcube.freeStorage (pcube, deleteIt1);
922 + fMat.putStorage (pflags, deleteIt2);
923 923 }
924 924
925 925
926 926 void VisImagingWeight::unPolChanWeight(Matrix<Float>& chanRowWt, const Cube<Float>& corrChanRowWt) {
927 927
928 928 //cout << "VIW::uPCW" << endl;
929 929
930 930 IPosition sh3=corrChanRowWt.shape();
931 931 IPosition sh2=sh3.getLast(2);
932 932 //cerr << "sh3" << sh3 << " sh2 " << sh2 << endl;

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut