Line data Source code
1 : //# VisImagingWeight.cc: imaging weight calculation for a give buffer
2 : //# Copyright (C) 2009-2018
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be adressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $Id$
28 : #include <msvis/MSVis/VisibilityIterator.h>
29 : #include <msvis/MSVis/VisBuffer.h>
30 : #include <msvis/MSVis/VisBuffer2.h>
31 : #include <msvis/MSVis/VisImagingWeight.h>
32 : #include <casacore/casa/Quanta/MVAngle.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/Matrix.h>
35 : #include <casacore/casa/Arrays/Vector.h>
36 : #include <casacore/lattices/Lattices/TempLattice.h>
37 : #include <casacore/images/Images/ImageInterface.h>
38 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
39 :
40 :
41 : using namespace casacore;
42 : namespace casa { //# NAMESPACE CASA - BEGIN
43 :
44 :
45 1526 : VisImagingWeight::VisImagingWeight() : wgtType_p("none"), doFilter_p(false), robust_p(0.0), rmode_p("norm"), noise_p(Quantity(0.0, "Jy")) , activeFieldIndex_p(-1){
46 :
47 1526 : }
48 :
49 141 : VisImagingWeight::VisImagingWeight(const String& type) : doFilter_p(false), robust_p(0.0), rmode_p("norm"), noise_p(Quantity(0.0, "Jy")) , activeFieldIndex_p(-1){
50 :
51 :
52 141 : wgtType_p=type;
53 141 : wgtType_p.downcase();
54 141 : if (wgtType_p != "natural" && wgtType_p != "radial"){
55 :
56 0 : throw(AipsError("Programmer error...wrong constructor used"));
57 : }
58 :
59 141 : }
60 :
61 :
62 0 : VisImagingWeight::VisImagingWeight(ROVisibilityIterator& vi, const String& rmode, const Quantity& noise,
63 : const Double robust, const Int nx, const Int ny,
64 : const Quantity& cellx, const Quantity& celly,
65 0 : const Int uBox, const Int vBox, const Bool multiField) : doFilter_p(false), robust_p(robust), rmode_p(rmode), noise_p(noise), activeFieldIndex_p(-1) {
66 :
67 0 : LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE));
68 :
69 0 : VisBufferAutoPtr vb (vi);
70 0 : wgtType_p="uniform";
71 : // Float uscale, vscale;
72 : //Int uorigin, vorigin;
73 0 : Vector<Double> deltas;
74 0 : uscale_p=(nx*cellx.get("rad").getValue());
75 0 : vscale_p=(ny*celly.get("rad").getValue());
76 0 : uorigin_p=nx/2;
77 0 : vorigin_p=ny/2;
78 0 : nx_p=nx;
79 0 : ny_p=ny;
80 : // Simply declare a big matrix
81 : //Matrix<Float> gwt(nx,ny);
82 0 : gwt_p.resize(1);
83 0 : multiFieldMap_p.clear();
84 0 : vi.originChunks();
85 0 : vi.origin();
86 :
87 : // Detect usable WEIGHT_SPECTRUM
88 0 : Bool doWtSp=(vb->weightSpectrum().nelements()>0);
89 :
90 0 : String mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId());
91 0 : multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(mapid, 0));
92 0 : gwt_p[0]=new TempLattice<Float>(IPosition(2, nx,ny), 0);
93 0 : gwt_p[0]->set(0.0);
94 0 : a_gwt_p.resize(nx_p, ny_p);
95 0 : a_gwt_p.set(0.0);
96 :
97 0 : Int fields=0;
98 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
99 0 : for (vi.origin();vi.more();vi++) {
100 :
101 0 : if(vb->newFieldId()){
102 0 : mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId());
103 0 : if(multiField){
104 0 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) ){
105 0 : fields+=1;
106 0 : gwt_p.resize(fields+1);
107 0 : gwt_p[fields]=new TempLattice<Float>(IPosition(2, nx_p,ny_p),0);
108 0 : gwt_p[fields]->set(0.0);
109 : }
110 : }
111 :
112 0 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) )
113 0 : multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(mapid, fields));
114 : }
115 : }
116 : }
117 :
118 : Float u, v;
119 0 : Vector<Double> sumwt(fields+1,0.0);
120 0 : f2_p.resize(fields+1);
121 0 : d2_p.resize(fields+1);
122 0 : Int fid=0;
123 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
124 0 : for (vi.origin();vi.more();vi++) {
125 0 : if(vb->newFieldId()){
126 0 : mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId());
127 :
128 :
129 0 : if(multiField){
130 0 : if(activeFieldIndex_p >=0)
131 0 : gwt_p[activeFieldIndex_p]->put(a_gwt_p);
132 : }
133 0 : activeFieldIndex_p=multiFieldMap_p.find(mapid) !=multiFieldMap_p.end( ) ? multiFieldMap_p[mapid] : -1;
134 0 : if(multiField && activeFieldIndex_p >=0)
135 0 : gwt_p[activeFieldIndex_p]->get(a_gwt_p);
136 :
137 : }
138 0 : auto mapvp = multiFieldMap_p.find(mapid);
139 0 : fid= mapvp != multiFieldMap_p.end( ) ? mapvp->second : -1;
140 0 : Int nRow=vb->nRow();
141 0 : Int nChan=vb->nChannel();
142 :
143 : // Extract weights correctly (use WEIGHT_SPECTRUM, if available)
144 0 : Matrix<Float> wtm;
145 : ////////////////
146 0 : doWtSp=(vb->weightSpectrum().nelements()>0);
147 : //////////////
148 0 : if (doWtSp)
149 : // WS available,
150 0 : unPolChanWeight(wtm,vb->weightSpectrum()); // collapse corr axis (parallel-hand average)
151 : else
152 : // WS UNavailable
153 0 : 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 0 : Int nChanWt=wtm.shape()(0);
157 :
158 0 : for (Int row=0; row<nRow; row++) {
159 0 : for (Int chn=0; chn<nChan; chn++) {
160 0 : if(!vb->flag()(chn,row)) {
161 0 : Float currwt=wtm(chn%nChanWt,row); // the weight for this chan,row
162 0 : Float f=vb->frequency()(chn)/C::c;
163 0 : u=vb->uvw()(row)(0)*f;
164 0 : v=vb->uvw()(row)(1)*f;
165 0 : Int ucell=Int(uscale_p*u+uorigin_p);
166 0 : Int vcell=Int(vscale_p*v+vorigin_p);
167 :
168 0 : if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) {
169 0 : for (Int iv=-vBox;iv<=vBox;iv++) {
170 0 : for (Int iu=-uBox;iu<=uBox;iu++) {
171 0 : a_gwt_p(ucell+iu,vcell+iv)+=currwt;
172 0 : sumwt[fid]+=currwt;
173 : }
174 : }
175 : }
176 0 : ucell=Int(-uscale_p*u+uorigin_p);
177 0 : vcell=Int(-vscale_p*v+vorigin_p);
178 0 : if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) {
179 0 : for (Int iv=-vBox;iv<=vBox;iv++) {
180 0 : for (Int iu=-uBox;iu<=uBox;iu++) {
181 0 : a_gwt_p(ucell+iu,vcell+iv)+=currwt;
182 0 : sumwt[fid]+=currwt;
183 : }
184 : }
185 : }
186 : }
187 : }
188 : }
189 0 : }
190 : }
191 : // make sure last active plane is saved
192 0 : gwt_p[activeFieldIndex_p]->put(a_gwt_p);
193 : // We use the approximation that all statistical weights are equal to
194 : // calculate the average summed weights (over visibilities, not bins!)
195 : // This is simply to try an ensure that the normalization of the robustness
196 : // parameter is similar to that of the ungridded case, but it doesn't have
197 : // to be exact, since any given case will require some experimentation.
198 :
199 : //Float f2, d2;
200 0 : for(fid=0; fid < Int(gwt_p.nelements()); ++fid){
201 0 : gwt_p[fid]->get(a_gwt_p);
202 0 : activeFieldIndex_p=fid;
203 0 : if (rmode=="norm") {
204 0 : os << "Normal robustness, robust = " << robust << LogIO::POST;
205 0 : Double sumlocwt = 0.;
206 0 : for(Int vgrid=0;vgrid<ny;vgrid++) {
207 0 : for(Int ugrid=0;ugrid<nx;ugrid++) {
208 0 : if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
209 : }
210 : }
211 0 : f2_p[fid]=0.0;
212 0 : if(sumlocwt !=0.0 && sumwt[fid] != 0.0)
213 0 : f2_p[fid] = square(5.0*pow(10.0,Double(-robust))) / (sumlocwt / sumwt[fid]);
214 0 : d2_p[fid] = 1.0;
215 :
216 : }
217 0 : else if (rmode=="abs") {
218 : os << "Absolute robustness, robust = " << robust << ", noise = "
219 0 : << noise.get("Jy").getValue() << "Jy" << LogIO::POST;
220 0 : f2_p[fid] = square(robust);
221 0 : d2_p[fid] = 2.0 * square(noise.get("Jy").getValue());
222 :
223 : }
224 : else {
225 0 : f2_p[fid] = 1.0;
226 0 : d2_p[fid] = 0.0;
227 : }
228 : }
229 :
230 0 : }
231 :
232 0 : VisImagingWeight::VisImagingWeight(ROVisibilityIterator& vi, Block<CountedPtr<TempLattice<Float> > >& grids, const String& rmode, const Quantity& noise,
233 : const Double robust, const Quantity& cellx, const Quantity& celly,
234 0 : const Bool multiField) : doFilter_p(false), robust_p(robust), rmode_p(rmode), noise_p(noise) {
235 :
236 0 : LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE));
237 :
238 0 : VisBufferAutoPtr vb (vi);
239 0 : wgtType_p="uniform";
240 0 : gwt_p.resize(grids.nelements(), true, false);
241 0 : for (uInt k =0; k < gwt_p.nelements(); ++k)
242 0 : gwt_p[k]=grids[k];
243 0 : nx_p=gwt_p[0]->shape()[0];
244 0 : ny_p=gwt_p[0]->shape()[1];
245 0 : uscale_p=(nx_p*cellx.get("rad").getValue());
246 0 : vscale_p=(ny_p*celly.get("rad").getValue());
247 0 : uorigin_p=nx_p/2;
248 0 : vorigin_p=ny_p/2;
249 :
250 0 : multiFieldMap_p.clear();
251 0 : Int fields=0;
252 0 : String mapid;
253 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
254 0 : for (vi.origin();vi.more();vi++) {
255 0 : if(vb->newFieldId()){
256 0 : mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId());
257 0 : if(multiField){
258 0 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) ){
259 0 : fields+=1;
260 :
261 : }
262 : }
263 0 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) )
264 0 : multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(mapid, fields));
265 : }
266 : }
267 : }
268 0 : f2_p.resize(gwt_p.nelements());
269 0 : d2_p.resize(gwt_p.nelements());
270 0 : for(uInt fid=0; fid < gwt_p.nelements(); ++fid){
271 0 : activeFieldIndex_p=fid;
272 0 : gwt_p[fid]->get(a_gwt_p);
273 0 : if (rmode_p=="norm") {
274 0 : os << "Normal robustness, robust = " << robust_p << LogIO::POST;
275 0 : Double sumlocwt = 0.;
276 0 : for(Int vgrid=0;vgrid<ny_p;vgrid++) {
277 0 : for(Int ugrid=0;ugrid<nx_p;ugrid++) {
278 0 : if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
279 : }
280 : }
281 0 : Double sumwt_fid=sum(a_gwt_p);
282 0 : f2_p[fid]=0.0;
283 0 : if(sumlocwt !=0.0 && sumwt_fid != 0)
284 0 : f2_p[fid] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumwt_fid);
285 0 : d2_p[fid] = 1.0;
286 :
287 : }
288 0 : else if (rmode=="abs") {
289 : os << "Absolute robustness, robust = " << robust_p << ", noise = "
290 0 : << noise_p.get("Jy").getValue() << "Jy" << LogIO::POST;
291 0 : f2_p[fid] = square(robust_p);
292 0 : d2_p[fid] = 2.0 * square(noise_p.get("Jy").getValue());
293 :
294 : }
295 : else {
296 0 : f2_p[fid] = 1.0;
297 0 : d2_p[fid] = 0.0;
298 : }
299 : }
300 0 : }
301 :
302 0 : VisImagingWeight::VisImagingWeight(vi::VisibilityIterator2& visIter,
303 : const String& rmode, const Quantity& noise,
304 : const Double robust, const Int nx, const Int ny,
305 : const Quantity& cellx, const Quantity& celly,
306 0 : const Int uBox, const Int vBox, const Bool multiField) : doFilter_p(false), robust_p(robust), rmode_p(rmode), noise_p(noise), activeFieldIndex_p(-1) {
307 :
308 0 : LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE));
309 :
310 :
311 :
312 0 : vi::VisBuffer2* vb=(visIter.getVisBuffer());
313 0 : wgtType_p="uniform";
314 : // Float uscale, vscale;
315 : //Int uorigin, vorigin;
316 0 : Vector<Double> deltas;
317 0 : uscale_p=(nx*cellx.get("rad").getValue());
318 0 : vscale_p=(ny*celly.get("rad").getValue());
319 0 : uorigin_p=nx/2;
320 0 : vorigin_p=ny/2;
321 0 : nx_p=nx;
322 0 : ny_p=ny;
323 : // Simply declare a big matrix
324 : //Matrix<Float> gwt(nx,ny);
325 0 : gwt_p.resize(1);
326 0 : multiFieldMap_p.clear();
327 0 : visIter.originChunks();
328 0 : visIter.origin();
329 0 : String mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]);
330 :
331 0 : multiFieldMap_p.insert( std::pair<casacore::String, casacore::Int>(mapid, 0) );
332 0 : gwt_p[0]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
333 0 : gwt_p[0]->set(0.0);
334 0 : a_gwt_p.resize(nx_p, ny_p);
335 0 : a_gwt_p.set(0.0);
336 :
337 :
338 : // Discover if weightSpectrum non-trivially available
339 0 : Bool doWtSp=visIter.weightSpectrumExists();
340 :
341 0 : Int fields=0;
342 0 : for (visIter.originChunks();visIter.moreChunks();visIter.nextChunk()) {
343 0 : for (visIter.origin();visIter.more();visIter.next()) {
344 0 : if(vb->isNewFieldId()){
345 0 : mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]);
346 0 : if(multiField){
347 0 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) ){
348 0 : fields+=1;
349 0 : gwt_p.resize(fields+1);
350 0 : gwt_p[fields]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
351 0 : gwt_p[fields]->set(0.0);
352 : }
353 : }
354 0 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) )
355 0 : multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(mapid, fields));
356 : }
357 : }
358 : }
359 :
360 : Float u, v;
361 0 : Vector<Double> sumwt(fields+1,0.0);
362 0 : f2_p.resize(fields+1);
363 0 : d2_p.resize(fields+1);
364 0 : Int fid=0;
365 0 : for (visIter.originChunks();visIter.moreChunks();visIter.nextChunk()) {
366 0 : for (visIter.origin();visIter.more();visIter.next()) {
367 0 : if(vb->isNewFieldId()){
368 0 : mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]);
369 :
370 0 : if(multiField){
371 0 : if(activeFieldIndex_p >=0)
372 0 : gwt_p[activeFieldIndex_p]->put(a_gwt_p);
373 : }
374 0 : activeFieldIndex_p= multiFieldMap_p.find(mapid) != multiFieldMap_p.end( ) ? multiFieldMap_p[mapid] : -1;
375 0 : if(multiField && activeFieldIndex_p >=0)
376 0 : gwt_p[activeFieldIndex_p]->get(a_gwt_p);
377 : }
378 :
379 0 : auto mapvp = multiFieldMap_p.find(mapid);
380 0 : fid= mapvp != multiFieldMap_p.end( ) ? mapvp->second : -1;
381 0 : Int nRow=vb->nRows();
382 0 : Int nChan=vb->nChannels();
383 :
384 : // Extract weights correctly (use WEIGHT_SPECTRUM, if available)
385 0 : Matrix<Float> wtm;
386 0 : Cube<Float> wtc;
387 0 : if (doWtSp)
388 : // WS available [nCorr,nChan,nRow]
389 0 : wtc.reference(vb->weightSpectrum());
390 : else {
391 0 : Int nCorr=vb->nCorrelations();
392 : // WS UNavailable weight()[nCorr,nRow] --> [nCorr,nChan,nRow]
393 0 : wtc.reference(vb->weight().reform(IPosition(3,nCorr,1,nRow))); // unchan'd weight as single-chan
394 : }
395 0 : unPolChanWeight(wtm,wtc); // Collapse on corr axis
396 :
397 : // Used for mod of chn index below
398 0 : Int nChanWt=wtm.shape()(0);
399 :
400 : //Oww !!! temporary implementation of old vb.flag just to see if things work
401 0 : Matrix<Bool> flag;
402 0 : cube2Matrix(vb->flagCube(), flag);
403 0 : for (Int row=0; row<nRow; row++) {
404 0 : for (Int chn=0; chn<nChan; chn++) {
405 :
406 0 : Float currwt=wtm(chn%nChanWt,row); // the weight for this chan,row
407 :
408 0 : if(!flag(chn,row)) {
409 0 : Float f=vb->getFrequency(row, chn)/C::c;
410 0 : u=vb->uvw()(0,row)*f;
411 0 : v=vb->uvw()(1,row)*f;
412 0 : Int ucell=Int(uscale_p*u+uorigin_p);
413 0 : Int vcell=Int(vscale_p*v+vorigin_p);
414 0 : if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) {
415 0 : for (Int iv=-vBox;iv<=vBox;iv++) {
416 0 : for (Int iu=-uBox;iu<=uBox;iu++) {
417 0 : a_gwt_p(ucell+iu,vcell+iv)+=currwt;
418 0 : sumwt[fid]+=currwt;
419 : }
420 : }
421 : }
422 0 : ucell=Int(-uscale_p*u+uorigin_p);
423 0 : vcell=Int(-vscale_p*v+vorigin_p);
424 0 : if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) {
425 0 : for (Int iv=-vBox;iv<=vBox;iv++) {
426 0 : for (Int iu=-uBox;iu<=uBox;iu++) {
427 0 : a_gwt_p(ucell+iu,vcell+iv)+=currwt;
428 0 : sumwt[fid]+=currwt;
429 : }
430 : }
431 : }
432 : }
433 : }
434 : }
435 0 : }
436 : }
437 :
438 : // make sure last active plane is saved
439 0 : gwt_p[activeFieldIndex_p]->put(a_gwt_p);
440 : // We use the approximation that all statistical weights are equal to
441 : // calculate the average summed weights (over visibilities, not bins!)
442 : // This is simply to try an ensure that the normalization of the robustness
443 : // parameter is similar to that of the ungridded case, but it doesn't have
444 : // to be exact, since any given case will require some experimentation.
445 :
446 : //Float f2, d2;
447 0 : for(fid=0; fid < Int(gwt_p.nelements()); ++fid){
448 0 : gwt_p[fid]->get(a_gwt_p);
449 0 : activeFieldIndex_p=fid;
450 0 : if (rmode=="norm") {
451 0 : os << "Normal robustness, robust = " << robust << LogIO::POST;
452 0 : Double sumlocwt = 0.;
453 0 : for(Int vgrid=0;vgrid<ny;vgrid++) {
454 0 : for(Int ugrid=0;ugrid<nx;ugrid++) {
455 0 : if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
456 : }
457 : }
458 0 : f2_p[fid]=0.0;
459 0 : if(sumlocwt != 0.0 && sumwt[fid] != 0.0)
460 0 : f2_p[fid] = square(5.0*pow(10.0,Double(-robust))) / (sumlocwt / sumwt[fid]);
461 :
462 0 : d2_p[fid] = 1.0;
463 :
464 : }
465 0 : else if (rmode=="abs") {
466 : os << "Absolute robustness, robust = " << robust << ", noise = "
467 0 : << noise.get("Jy").getValue() << "Jy" << LogIO::POST;
468 0 : f2_p[fid] = square(robust);
469 0 : d2_p[fid] = 2.0 * square(noise.get("Jy").getValue());
470 :
471 : }
472 : else {
473 0 : f2_p[fid] = 1.0;
474 0 : d2_p[fid] = 0.0;
475 : }
476 : }
477 0 : }
478 :
479 :
480 0 : VisImagingWeight::VisImagingWeight(ImageInterface<Float>& im) : robust_p(0.0), rmode_p(""), noise_p(Quantity(0.0, "Jy")) {
481 :
482 0 : LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE));
483 :
484 0 : doFilter_p=False;
485 :
486 0 : wgtType_p="uniform";
487 0 : nx_p=im.shape()(0);
488 0 : ny_p=im.shape()(1);
489 0 : DirectionCoordinate dc=im.coordinates().directionCoordinate(0);
490 0 : dc.setWorldAxisUnits(Vector<String>(2, "rad"));
491 0 : Double cellx=dc.increment()(0);
492 0 : Double celly=dc.increment()(1);
493 0 : uscale_p=nx_p*cellx;
494 0 : vscale_p=ny_p*celly;
495 0 : uorigin_p=nx_p/2;
496 0 : vorigin_p=ny_p/2;
497 : //Now to recover from image density and other parameters
498 0 : Int nplanes=1;
499 0 : if(im.shape().nelements()==5)
500 0 : nplanes=im.shape()[4];
501 0 : gwt_p.resize(nplanes, True, False);
502 0 : if(im.shape().nelements()==5){
503 0 : IPosition blc(Vector<Int>(5,0));
504 :
505 0 : for (Int fid=0;fid<nplanes;fid++)
506 : {
507 0 : gwt_p[fid]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
508 0 : Array<Float> lala;
509 0 : blc[4]=fid;
510 0 : im.getSlice(lala, blc, IPosition(5, nx_p, ny_p,1,1,1), True);
511 0 : gwt_p[fid]->put( lala.reform(IPosition(2, nx_p, ny_p)));
512 :
513 0 : }
514 0 : }
515 : else{
516 0 : Array<Float> lala;
517 0 : im.get(lala, True);
518 0 : gwt_p[0]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
519 0 : gwt_p[0]->put( lala.reform(IPosition(2, nx_p, ny_p)));
520 :
521 0 : }
522 0 : const TableRecord& rec=im.miscInfo();
523 0 : if(rec.isDefined("d2")){
524 0 : d2_p.resize();
525 0 : rec.get("d2", d2_p);
526 0 : f2_p.resize();
527 0 : rec.get("f2", f2_p);
528 0 : multiFieldMap_p.clear();
529 : Int mapsize;
530 0 : rec.get("multimapsize", mapsize);
531 0 : for(Int k=0; k < mapsize; ++k){
532 0 : String key;
533 : Int val;
534 0 : rec.get("key"+String::toString(k), key);
535 0 : rec.get("val"+String::toString(k), val);
536 0 : multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(key, val));
537 0 : }
538 :
539 :
540 : }
541 0 : activeFieldIndex_p=0;
542 0 : gwt_p[0]->get(a_gwt_p);
543 0 : if(rec.isDefined("dofilter")){
544 0 : rec.get("dofilter", doFilter_p);
545 0 : rec.get("rbmaj", rbmaj_p);
546 0 : rec.get("rbmin", rbmin_p);
547 0 : rec.get("cospa", cospa_p);
548 0 : rec.get("sinpa", sinpa_p);
549 : }
550 :
551 :
552 0 : }
553 :
554 :
555 1665 : VisImagingWeight::~VisImagingWeight(){
556 : /*for (uInt fid=0; fid < gwt_p.nelements(); ++fid){
557 : gwt_p[fid].resize();
558 : }*/
559 :
560 1665 : }
561 :
562 0 : Vector<Int> VisImagingWeight::shapeOfdensityGrid(){
563 0 : Vector<Int> retval(3, 0);
564 0 : retval(2)=gwt_p.nelements();
565 0 : if(retval(2) > 0){
566 0 : retval[0]=gwt_p[0]->shape()(0);
567 0 : retval[1]=gwt_p[0]->shape()(1);
568 :
569 : }
570 :
571 0 : return retval;
572 0 : }
573 0 : void VisImagingWeight::toImageInterface(casacore::ImageInterface<casacore::Float>& im){
574 :
575 0 : if( wgtType_p != "uniform")
576 0 : throw(AipsError("cannot save weight density for non-Briggs' weighting schemes"));
577 :
578 0 : IPosition where= IPosition(im.shape().nelements(),0);
579 0 : Int lastAx=where.nelements()-1;
580 0 : for (uInt fid=0;fid<gwt_p.nelements(); ++fid){
581 0 : activeFieldIndex_p=fid;
582 0 : gwt_p[fid]->get(a_gwt_p);
583 0 : where[lastAx]=fid;
584 0 : im.putSlice(a_gwt_p, where);
585 :
586 : }
587 0 : Record rec(im.miscInfo());
588 0 : rec.define("d2", d2_p);
589 0 : rec.define("f2", f2_p);
590 0 : rec.define("numfield", Int(multiFieldMap_p.size()));
591 0 : uInt keycount=0;
592 0 : for (auto iter = multiFieldMap_p.begin( ); iter != multiFieldMap_p.end( ); ++iter, ++keycount){
593 0 : rec.define("key"+String::toString(keycount), iter->first);
594 0 : rec.define("val"+String::toString(keycount), iter->second);
595 : }
596 0 : rec.define("multimapsize",keycount);
597 0 : if(doFilter_p){
598 0 : rec.define("dofilter", doFilter_p);
599 0 : rec.define("rbmaj", rbmaj_p);
600 0 : rec.define("rbmin", rbmin_p);
601 0 : rec.define("cospa", cospa_p);
602 0 : rec.define("sinpa", sinpa_p);
603 : }
604 :
605 :
606 0 : im.setMiscInfo(rec);
607 :
608 0 : }
609 0 : void VisImagingWeight::setFilter(const String& type, const Quantity& bmaj,
610 : const Quantity& bmin, const Quantity& bpa)
611 : {
612 :
613 0 : LogIO os(LogOrigin("VisImagingWeight", "setFilter()", WHERE));
614 :
615 : // Steps : (1) Convert the uvtaper information into "sigma" for the uv-domain Gaussian.
616 : // (2) Rotate u,v per cell in the weightdensity grid onto an axis where the position angle is zero. The u,v becomes ru,rv.
617 : // (3) Evaluate the Gaussian as e^{-1/2 (ru^2/sigma_{maj}^2 + rv^2/sigma_{min}^2)} where sigma is separate for the maj and min axes.
618 : // (4) Multiply the weightdensity grid (each cell) with the evaluated Gaussian above.
619 :
620 : // There are two ways this input may be supplied : In the image domain or the UV domain.
621 : // For both options, the code below calculates xx = 1 / ( 2 sigma_{maj}^2) and yy = 1 / (2 sigma_{min}^2) for the major and minor axes.
622 : // The Gaussian is then evaluated as exp( - xx * ru^2 - yy * rv^2 ) in the ::filter() method. In this code, xx is called rbmaj_p and rbmin_p.
623 :
624 0 : if (type=="gaussian") {
625 :
626 : // uvtaper input is supplied in the uv domain as the HWHM of a Gaussian. This is a 'uvdistance' or 'baseline length' interpretation
627 : // The FWHM_uv (in wavelengths) = beam_lambda * 2 to take it from HWHM to FWHM.
628 : // sigma_uv = FWHM_uv / (2 sqrt(2 log2))
629 : // xx = 1 / ( 2 sigma_uv^2) = (log2) / (beam_lambda)^2
630 :
631 0 : Bool lambdafilt=false;
632 :
633 0 : if( bmaj.getUnit().contains("lambda"))
634 0 : lambdafilt=true;
635 0 : if(lambdafilt){
636 : os << "Filtering for Gaussian of shape from read: "
637 0 : << bmaj.get("klambda").getValue() << " by "
638 0 : << bmin.get("klambda").getValue() << " (klambda) at p.a. "
639 0 : << bpa.get("deg").getValue() << " (degrees)" << LogIO::POST;
640 :
641 0 : rbmaj_p=log(2.0)/square(bmaj.get("lambda").getValue());
642 0 : rbmin_p=log(2.0)/square(bmin.get("lambda").getValue());
643 :
644 : }
645 : else{
646 :
647 : // uvtaper input is supplied in the image domain, as the FWHM of a Gaussian. This is an 'convolving' angular resolution interpretation
648 : // This FWHM_lm (in radians) must be converted to a "Sigma" of the Gaussian, and then taken to the UV-domain.
649 : // FWHM_lm = beam_radians
650 : // sigma_lm = FWHM_lm / (2 sqrt(2 log2))
651 : // sigma_uv = 1 / ( 2 pi sigma_lm )
652 : // xx = 1 / ( 2 sigma_uv^2) = ( (pi^2)/(4 log2) ) * (beam_radians)^2
653 :
654 : os << "Filtering for Gaussian of shape: "
655 0 : << bmaj.get("arcsec").getValue() << " by "
656 0 : << bmin.get("arcsec").getValue() << " (arcsec) at p.a. "
657 0 : << bpa.get("deg").getValue() << " (degrees)" << LogIO::POST;
658 :
659 : // Convert to values that we can use
660 0 : Double fact = M_PI*M_PI/(4.0*log(2.0));
661 0 : rbmaj_p = fact*square(bmaj.get("rad").getValue());
662 0 : rbmin_p = fact*square(bmin.get("rad").getValue());
663 : }
664 :
665 0 : Double rbpa = MVAngle(bpa).get("rad").getValue();
666 0 : cospa_p = sin(rbpa);
667 0 : sinpa_p = cos(rbpa);
668 :
669 : os << "Filtering for Gaussian of shape after convention: maj"
670 : << rbmaj_p << " min "
671 : << rbmin_p<< " pa "
672 0 : << rbpa << " " << LogIO::POST;
673 :
674 0 : doFilter_p=true;
675 : }
676 : else {
677 0 : os << "Unknown filtering " << type << LogIO::EXCEPTION;
678 : }
679 :
680 0 : }
681 :
682 :
683 6520 : Bool VisImagingWeight::doFilter() const{
684 :
685 6520 : return doFilter_p;
686 : }
687 :
688 :
689 0 : void VisImagingWeight::filter(Matrix<Float>& imWeight, const Matrix<Bool>& flag,
690 : const Matrix<Double>& uvw,
691 : const Vector<Double>& frequency, const Matrix<Float>& weight) const{
692 :
693 : // Expecting weight[nchan,nrow], where nchan=1 or nChan (data)
694 : // If nchan=1 (i.e., WEIGHT_SPECTRUM unavailable), then the
695 : // following will be handle the channel degeneracy correctly, using %.
696 :
697 :
698 0 : Int nRow=imWeight.shape()(1);
699 0 : Int nChan=imWeight.shape()(0);
700 0 : Int nChanWt=weight.shape()(0);
701 0 : for (Int row=0; row<nRow; row++) {
702 0 : for (Int chn=0; chn<nChan; chn++) {
703 0 : Double invLambdaC=frequency(chn)/C::c;
704 0 : Double u = uvw(0,row);
705 0 : Double v = uvw(1,row);
706 0 : if(!flag(chn,row) && (weight(chn%nChanWt,row)>0.0) ) {
707 : // Rotate the u,v values of each cell, so that 'ru' and 'rv' are aligned with the Major and Minor axie of the uvtaper Gaussian.
708 : // If the uvtaper Gaussian has positionangle=0, then ru=u, rv=v.
709 : // If the uvtaper Gaussian has positionangle=90, then ru=v, rv= -u
710 0 : Double ru = invLambdaC*( cospa_p * u + sinpa_p * v);
711 0 : Double rv = invLambdaC*(- sinpa_p * u + cospa_p * v);
712 0 : Double filter = exp(-rbmaj_p*square(ru) - rbmin_p*square(rv));
713 0 : imWeight(chn,row)*=filter;
714 : }
715 : else {
716 0 : imWeight(chn,row)=0.0;
717 : }
718 : }
719 : }
720 :
721 :
722 0 : }
723 :
724 :
725 0 : void VisImagingWeight::weightUniform(Matrix<Float>& imWeight, const Matrix<Bool>& flag, const Matrix<Double>& uvw,
726 : const Vector<Double>& frequency,
727 : const Matrix<Float>& weight, const Int msId, const Int fieldId) const{
728 :
729 :
730 : //cout << " WEIG " << nx_p << " " << ny_p << " " << gwt_p[0].shape() << endl;
731 : //cout << "f2 " << f2_p[0] << " d2 " << d2_p[0] << " uscale " << uscale_p << " vscal " << vscale_p << " origs " << uorigin_p << " " << vorigin_p << endl;
732 0 : String mapid=String::toString(msId)+String("_")+String::toString(fieldId);
733 : //cout << "min max gwt " << min(gwt_p[0]) << " " << max(gwt_p[0]) << " mapid " << mapid <<endl;
734 :
735 0 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) )
736 0 : throw(AipsError("Imaging weight calculation is requested for a data that was not selected"));
737 :
738 0 : auto mapvp = multiFieldMap_p.find(mapid);
739 0 : Int fid = mapvp != multiFieldMap_p.end( ) ? mapvp->second : -1;
740 : //Int ndrop=0;
741 0 : if(activeFieldIndex_p != fid){
742 0 : a_gwt_p=gwt_p[fid]->get();
743 0 : activeFieldIndex_p=fid;
744 : }
745 0 : Double sumwt=0.0;
746 0 : Int nRow=imWeight.shape()(1);
747 0 : Int nChannel=imWeight.shape()(0);
748 0 : Int nChanWt=weight.shape()(0);
749 : Float u, v;
750 0 : for (Int row=0; row<nRow; row++) {
751 0 : for (Int chn=0; chn<nChannel; chn++) {
752 0 : if (!flag(chn,row)) {
753 0 : Float f=frequency(chn)/C::c;
754 0 : u=uvw(0, row)*f;
755 0 : v=uvw(1, row)*f;
756 0 : Int ucell=Int(uscale_p*u+uorigin_p);
757 0 : Int vcell=Int(vscale_p*v+vorigin_p);
758 0 : imWeight(chn,row)=weight(chn%nChanWt,row);
759 0 : if((ucell>0)&&(ucell<nx_p)&&(vcell>0)&&(vcell<ny_p) &&a_gwt_p(ucell,vcell)>0.0) {
760 0 : imWeight(chn,row)/=a_gwt_p(ucell,vcell)*f2_p[fid]+d2_p[fid];
761 0 : sumwt+=imWeight(chn,row);
762 :
763 : }
764 : else {
765 0 : imWeight(chn,row)=0.0;
766 : //ndrop++;
767 : }
768 : }
769 : else{
770 0 : imWeight(chn,row)=0.0;
771 : }
772 : }
773 : }
774 :
775 0 : }
776 :
777 : /* unused version?
778 : void VisImagingWeight::weightNatural(Matrix<Float>& imagingWeight, const Matrix<Bool>& flag,
779 : const Matrix<Float>& weight) const{
780 :
781 :
782 :
783 : Int nRow=imagingWeight.shape()(1);
784 : Vector<Float> wgtRow(nRow);
785 :
786 : for (Int row=0; row<nRow; row++) {
787 : wgtRow(row)=max(weight.column(row));
788 : }
789 : weightNatural(imagingWeight, flag, wgtRow);
790 :
791 : }
792 : */
793 6320 : void VisImagingWeight::weightNatural(Matrix<Float>& imagingWeight, const Matrix<Bool>& flag,
794 : const Matrix<Float>& weight) const{
795 :
796 6320 : Double sumwt=0.0;
797 : //cerr << "shape of weight " << weight.shape() << " max " << max (weight.column(0) ) << " max " << min(weight.column(0)) << " " << weight.column(0).shape() << endl;
798 6320 : Int nRow=imagingWeight.shape()(1);
799 6320 : Int nChan=imagingWeight.shape()(0);
800 6320 : Int nChanWt=weight.shape()(0);
801 810920 : for (Int row=0; row<nRow; row++) {
802 8559000 : for (Int chn=0; chn<nChan; chn++) {
803 7754400 : if( !flag(chn,row) ) {
804 7754400 : imagingWeight(chn,row)=weight(chn%nChanWt,row);
805 7754400 : sumwt+=imagingWeight(chn,row);
806 : }
807 : else {
808 0 : imagingWeight(chn,row)=0.0;
809 : }
810 : }
811 : }
812 :
813 :
814 6320 : }
815 :
816 :
817 0 : void VisImagingWeight::weightRadial(Matrix<Float>& imagingWeight,
818 : const Matrix<Bool>& flag,
819 : const Matrix<Double>& uvw,
820 : const Vector<Double>& frequency,
821 : const Matrix<Float>& weight) const{
822 :
823 0 : Double sumwt=0.0;
824 0 : Int nRow=imagingWeight.shape()(1);
825 0 : Int nChan=imagingWeight.shape()(0);
826 0 : Int nChanWt=weight.shape()(0);
827 :
828 0 : for (Int row=0; row<nRow; row++) {
829 0 : for (Int chn=0; chn< nChan; chn++) {
830 0 : Float f=frequency(chn)/C::c;
831 0 : if( !flag(chn,row) ) {
832 0 : imagingWeight(chn,row)=
833 0 : f*sqrt(square(uvw(0, row))+square(uvw(1, row)))
834 0 : *weight(chn%nChanWt,row);
835 0 : sumwt+=imagingWeight(chn,row);
836 : }
837 : else {
838 0 : imagingWeight(chn,row)=0.0;
839 : }
840 : }
841 : }
842 :
843 :
844 0 : }
845 :
846 529 : VisImagingWeight& VisImagingWeight::operator=(const VisImagingWeight& other){
847 529 : if(this != &other){
848 : // gwt_p=other.gwt_p;
849 :
850 529 : gwt_p.resize(other.gwt_p.nelements(), true, false);
851 529 : for (uInt k=0; k < gwt_p.nelements(); ++k){
852 : //gwt_p[k].resize();
853 0 : gwt_p[k]=other.gwt_p[k];
854 : }
855 :
856 :
857 529 : wgtType_p=other.wgtType_p;
858 529 : uscale_p=other.uscale_p;
859 529 : vscale_p=other.vscale_p;
860 529 : f2_p.resize();
861 529 : d2_p.resize();
862 529 : f2_p=other.f2_p;
863 529 : d2_p=other.d2_p;
864 529 : uorigin_p=other.uorigin_p;
865 529 : vorigin_p=other.vorigin_p;
866 529 : nx_p=other.nx_p;
867 529 : ny_p=other.ny_p;
868 529 : doFilter_p=other.doFilter_p;
869 529 : cospa_p=other.cospa_p;
870 529 : sinpa_p=other.sinpa_p;
871 529 : rbmaj_p=other.rbmaj_p;
872 529 : rbmin_p=other.rbmin_p;
873 529 : robust_p=other.robust_p;
874 529 : rmode_p=other.rmode_p;
875 529 : multiFieldMap_p=other.multiFieldMap_p;
876 : }
877 529 : return *this;
878 : }
879 :
880 19266 : String VisImagingWeight::getType() const{
881 :
882 19266 : return wgtType_p;
883 :
884 : }
885 0 : Bool VisImagingWeight::getWeightDensity (Block<Matrix<Float> >& density){
886 0 : if(wgtType_p != "uniform"){
887 0 : density.resize(0, true, false);
888 0 : return false;
889 : }
890 0 : density.resize(gwt_p.nelements(), true, false);
891 0 : for (uInt k=0; k < gwt_p.nelements(); ++k){
892 0 : density[k].resize(gwt_p[k]->shape());
893 0 : density[k]=gwt_p[k]->get();
894 : }
895 0 : return true;
896 : }
897 0 : void VisImagingWeight::setWeightDensity(const Block<Matrix<Float> >& density){
898 0 : if(wgtType_p=="uniform"){
899 0 : gwt_p.resize(density.nelements(), true, false);
900 0 : for (uInt k=0; k < gwt_p.nelements(); ++k){
901 : //gwt_p[k].resize();
902 0 : gwt_p[k]=new TempLattice<Float>(IPosition(2, density[k].shape()[0], density[k].shape()[1]),0);
903 0 : gwt_p[k]->put(density[k]);
904 : }
905 : //break any old reference
906 0 : f2_p.resize();
907 0 : d2_p.resize();
908 0 : f2_p.resize(gwt_p.nelements());
909 0 : d2_p.resize(gwt_p.nelements());
910 0 : f2_p.set(1.0);
911 0 : d2_p.set(0.0);
912 : //Float f2, d2;
913 0 : for(uInt fid=0; fid < gwt_p.nelements(); ++fid){
914 0 : activeFieldIndex_p=fid;
915 0 : a_gwt_p=gwt_p[fid]->get();
916 0 : if (rmode_p=="norm") {
917 0 : Double sumlocwt = 0.;
918 0 : for(Int vgrid=0;vgrid<a_gwt_p.shape()(1);vgrid++) {
919 0 : for(Int ugrid=0;ugrid<a_gwt_p.shape()(0);ugrid++) {
920 0 : if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
921 : }
922 : }
923 0 : Double sumwt_fid=sum(a_gwt_p);
924 0 : f2_p[fid]=0.0;
925 0 : if(sumlocwt != 0.0 && sumwt_fid != 0.0)
926 0 : f2_p[fid] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumwt_fid);
927 0 : d2_p[fid] = 1.0;
928 : }
929 0 : else if (rmode_p=="abs") {
930 0 : f2_p[fid] = square(robust_p);
931 0 : d2_p[fid] = 2.0 * square(noise_p.get("Jy").getValue());
932 :
933 : }
934 : else {
935 0 : f2_p[fid] = 1.0;
936 0 : d2_p[fid] = 0.0;
937 : }
938 : }
939 :
940 : }
941 :
942 :
943 0 : }
944 0 : void VisImagingWeight::cube2Matrix(const Cube<Bool>& fcube, Matrix<Bool>& fMat)
945 : {
946 0 : fMat.resize(fcube.shape()[1], fcube.shape()[2]);
947 : Bool deleteIt1;
948 : Bool deleteIt2;
949 0 : const Bool * pcube = fcube.getStorage (deleteIt1);
950 0 : Bool * pflags = fMat.getStorage (deleteIt2);
951 0 : for (uInt row = 0; row < fcube.shape()[2]; row++) {
952 0 : for (Int chn = 0; chn < fcube.shape()[1]; chn++) {
953 0 : *pflags = *pcube++;
954 0 : for (Int pol = 1; pol < fcube.shape()[0]; pol++, pcube++) {
955 0 : *pflags = *pcube ? *pcube : *pflags;
956 : }
957 0 : pflags++;
958 : }
959 : }
960 0 : fcube.freeStorage (pcube, deleteIt1);
961 0 : fMat.putStorage (pflags, deleteIt2);
962 0 : }
963 :
964 :
965 6120 : void VisImagingWeight::unPolChanWeight(Matrix<Float>& chanRowWt, const Cube<Float>& corrChanRowWt) {
966 :
967 : //cout << "VIW::uPCW" << endl;
968 :
969 6120 : IPosition sh3=corrChanRowWt.shape();
970 6120 : IPosition sh2=sh3.getLast(2);
971 : //cerr << "sh3" << sh3 << " sh2 " << sh2 << endl;
972 6120 : chanRowWt.resize(sh2);
973 : //cout << "assign..." << endl;
974 6120 : chanRowWt=corrChanRowWt(Slice(0,1,1),Slice(),Slice()).reform(sh2);
975 6120 : const Int& nPol=sh3[0];
976 6120 : if (nPol>1) {
977 6120 : chanRowWt += corrChanRowWt(Slice(nPol-1,1,1),Slice(),Slice()).reform(sh2);
978 6120 : chanRowWt /= 2.0f;
979 : }
980 :
981 : // cout << "VIW::uPCW" << endl;
982 :
983 6120 : }
984 :
985 :
986 : }//# NAMESPACE CASA - END
987 :
988 :
|
