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 12444 : 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 12444 : }
48 :
49 4277 : 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 4277 : wgtType_p=type;
53 4277 : wgtType_p.downcase();
54 4277 : if (wgtType_p != "natural" && wgtType_p != "radial"){
55 :
56 0 : throw(AipsError("Programmer error...wrong constructor used"));
57 : }
58 :
59 4277 : }
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] = square(5.0*pow(10.0,Double(-robust))) / (sumlocwt / sumwt[fid]);
212 0 : d2_p[fid] = 1.0;
213 :
214 : }
215 0 : else if (rmode=="abs") {
216 : os << "Absolute robustness, robust = " << robust << ", noise = "
217 0 : << noise.get("Jy").getValue() << "Jy" << LogIO::POST;
218 0 : f2_p[fid] = square(robust);
219 0 : d2_p[fid] = 2.0 * square(noise.get("Jy").getValue());
220 :
221 : }
222 : else {
223 0 : f2_p[fid] = 1.0;
224 0 : d2_p[fid] = 0.0;
225 : }
226 : }
227 :
228 0 : }
229 :
230 0 : VisImagingWeight::VisImagingWeight(ROVisibilityIterator& vi, Block<CountedPtr<TempLattice<Float> > >& grids, const String& rmode, const Quantity& noise,
231 : const Double robust, const Quantity& cellx, const Quantity& celly,
232 0 : const Bool multiField) : doFilter_p(false), robust_p(robust), rmode_p(rmode), noise_p(noise) {
233 :
234 0 : LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE));
235 :
236 0 : VisBufferAutoPtr vb (vi);
237 0 : wgtType_p="uniform";
238 0 : gwt_p.resize(grids.nelements(), true, false);
239 0 : for (uInt k =0; k < gwt_p.nelements(); ++k)
240 0 : gwt_p[k]=grids[k];
241 0 : nx_p=gwt_p[0]->shape()[0];
242 0 : ny_p=gwt_p[0]->shape()[1];
243 0 : uscale_p=(nx_p*cellx.get("rad").getValue());
244 0 : vscale_p=(ny_p*celly.get("rad").getValue());
245 0 : uorigin_p=nx_p/2;
246 0 : vorigin_p=ny_p/2;
247 :
248 0 : multiFieldMap_p.clear();
249 0 : Int fields=0;
250 0 : String mapid;
251 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
252 0 : for (vi.origin();vi.more();vi++) {
253 0 : if(vb->newFieldId()){
254 0 : mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId());
255 0 : if(multiField){
256 0 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) ){
257 0 : fields+=1;
258 :
259 : }
260 : }
261 0 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) )
262 0 : multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(mapid, fields));
263 : }
264 : }
265 : }
266 0 : f2_p.resize(gwt_p.nelements());
267 0 : d2_p.resize(gwt_p.nelements());
268 0 : for(uInt fid=0; fid < gwt_p.nelements(); ++fid){
269 0 : activeFieldIndex_p=fid;
270 0 : gwt_p[fid]->get(a_gwt_p);
271 0 : if (rmode_p=="norm") {
272 0 : os << "Normal robustness, robust = " << robust_p << LogIO::POST;
273 0 : Double sumlocwt = 0.;
274 0 : for(Int vgrid=0;vgrid<ny_p;vgrid++) {
275 0 : for(Int ugrid=0;ugrid<nx_p;ugrid++) {
276 0 : if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
277 : }
278 : }
279 0 : Double sumwt_fid=sum(a_gwt_p);
280 0 : f2_p[fid] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumwt_fid);
281 0 : d2_p[fid] = 1.0;
282 :
283 : }
284 0 : else if (rmode=="abs") {
285 : os << "Absolute robustness, robust = " << robust_p << ", noise = "
286 0 : << noise_p.get("Jy").getValue() << "Jy" << LogIO::POST;
287 0 : f2_p[fid] = square(robust_p);
288 0 : d2_p[fid] = 2.0 * square(noise_p.get("Jy").getValue());
289 :
290 : }
291 : else {
292 0 : f2_p[fid] = 1.0;
293 0 : d2_p[fid] = 0.0;
294 : }
295 : }
296 0 : }
297 :
298 18 : VisImagingWeight::VisImagingWeight(vi::VisibilityIterator2& visIter,
299 : const String& rmode, const Quantity& noise,
300 : const Double robust, const Int nx, const Int ny,
301 : const Quantity& cellx, const Quantity& celly,
302 18 : 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 :
304 36 : LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE));
305 :
306 :
307 :
308 18 : vi::VisBuffer2* vb=(visIter.getVisBuffer());
309 18 : wgtType_p="uniform";
310 : // Float uscale, vscale;
311 : //Int uorigin, vorigin;
312 18 : Vector<Double> deltas;
313 18 : uscale_p=(nx*cellx.get("rad").getValue());
314 18 : vscale_p=(ny*celly.get("rad").getValue());
315 18 : uorigin_p=nx/2;
316 18 : vorigin_p=ny/2;
317 18 : nx_p=nx;
318 18 : ny_p=ny;
319 : // Simply declare a big matrix
320 : //Matrix<Float> gwt(nx,ny);
321 18 : gwt_p.resize(1);
322 18 : multiFieldMap_p.clear();
323 18 : visIter.originChunks();
324 18 : visIter.origin();
325 36 : String mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]);
326 :
327 18 : multiFieldMap_p.insert( std::pair<casacore::String, casacore::Int>(mapid, 0) );
328 18 : gwt_p[0]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
329 18 : gwt_p[0]->set(0.0);
330 18 : a_gwt_p.resize(nx_p, ny_p);
331 18 : a_gwt_p.set(0.0);
332 :
333 :
334 : // Discover if weightSpectrum non-trivially available
335 18 : Bool doWtSp=visIter.weightSpectrumExists();
336 :
337 18 : Int fields=0;
338 51 : for (visIter.originChunks();visIter.moreChunks();visIter.nextChunk()) {
339 5793 : for (visIter.origin();visIter.more();visIter.next()) {
340 5760 : if(vb->isNewFieldId()){
341 5712 : mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]);
342 5712 : if(multiField){
343 1152 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) ){
344 9 : fields+=1;
345 9 : gwt_p.resize(fields+1);
346 9 : gwt_p[fields]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
347 9 : gwt_p[fields]->set(0.0);
348 : }
349 : }
350 5712 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) )
351 11 : multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(mapid, fields));
352 : }
353 : }
354 : }
355 :
356 : Float u, v;
357 18 : Vector<Double> sumwt(fields+1,0.0);
358 18 : f2_p.resize(fields+1);
359 18 : d2_p.resize(fields+1);
360 18 : Int fid=0;
361 51 : for (visIter.originChunks();visIter.moreChunks();visIter.nextChunk()) {
362 5793 : for (visIter.origin();visIter.more();visIter.next()) {
363 5760 : if(vb->isNewFieldId()){
364 5712 : mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]);
365 :
366 5712 : if(multiField){
367 1152 : if(activeFieldIndex_p >=0)
368 1149 : gwt_p[activeFieldIndex_p]->put(a_gwt_p);
369 : }
370 5712 : activeFieldIndex_p= multiFieldMap_p.find(mapid) != multiFieldMap_p.end( ) ? multiFieldMap_p[mapid] : -1;
371 5712 : if(multiField && activeFieldIndex_p >=0)
372 1152 : gwt_p[activeFieldIndex_p]->get(a_gwt_p);
373 : }
374 :
375 5760 : auto mapvp = multiFieldMap_p.find(mapid);
376 5760 : fid= mapvp != multiFieldMap_p.end( ) ? mapvp->second : -1;
377 5760 : Int nRow=vb->nRows();
378 5760 : Int nChan=vb->nChannels();
379 :
380 : // Extract weights correctly (use WEIGHT_SPECTRUM, if available)
381 5760 : Matrix<Float> wtm;
382 5760 : Cube<Float> wtc;
383 5760 : if (doWtSp)
384 : // WS available [nCorr,nChan,nRow]
385 0 : wtc.reference(vb->weightSpectrum());
386 : else {
387 5760 : Int nCorr=vb->nCorrelations();
388 : // WS UNavailable weight()[nCorr,nRow] --> [nCorr,nChan,nRow]
389 5760 : wtc.reference(vb->weight().reform(IPosition(3,nCorr,1,nRow))); // unchan'd weight as single-chan
390 : }
391 5760 : unPolChanWeight(wtm,wtc); // Collapse on corr axis
392 :
393 : // Used for mod of chn index below
394 5760 : Int nChanWt=wtm.shape()(0);
395 :
396 : //Oww !!! temporary implementation of old vb.flag just to see if things work
397 5760 : Matrix<Bool> flag;
398 5760 : cube2Matrix(vb->flagCube(), flag);
399 1958516 : for (Int row=0; row<nRow; row++) {
400 11326992 : for (Int chn=0; chn<nChan; chn++) {
401 :
402 9374236 : Float currwt=wtm(chn%nChanWt,row); // the weight for this chan,row
403 :
404 9374236 : if(!flag(chn,row)) {
405 8972041 : Float f=vb->getFrequency(row, chn)/C::c;
406 8972041 : u=vb->uvw()(0,row)*f;
407 8972041 : v=vb->uvw()(1,row)*f;
408 8972041 : Int ucell=Int(uscale_p*u+uorigin_p);
409 8972041 : Int vcell=Int(vscale_p*v+vorigin_p);
410 8972041 : if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) {
411 19874428 : for (Int iv=-vBox;iv<=vBox;iv++) {
412 35320156 : for (Int iu=-uBox;iu<=uBox;iu++) {
413 24417584 : a_gwt_p(ucell+iu,vcell+iv)+=currwt;
414 24417584 : sumwt[fid]+=currwt;
415 : }
416 : }
417 : }
418 8972041 : ucell=Int(-uscale_p*u+uorigin_p);
419 8972041 : vcell=Int(-vscale_p*v+vorigin_p);
420 8972041 : if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) {
421 19874434 : for (Int iv=-vBox;iv<=vBox;iv++) {
422 35320162 : for (Int iu=-uBox;iu<=uBox;iu++) {
423 24417587 : a_gwt_p(ucell+iu,vcell+iv)+=currwt;
424 24417587 : sumwt[fid]+=currwt;
425 : }
426 : }
427 : }
428 : }
429 : }
430 : }
431 5760 : }
432 : }
433 :
434 : // make sure last active plane is saved
435 18 : gwt_p[activeFieldIndex_p]->put(a_gwt_p);
436 : // We use the approximation that all statistical weights are equal to
437 : // calculate the average summed weights (over visibilities, not bins!)
438 : // This is simply to try an ensure that the normalization of the robustness
439 : // parameter is similar to that of the ungridded case, but it doesn't have
440 : // to be exact, since any given case will require some experimentation.
441 :
442 : //Float f2, d2;
443 45 : for(fid=0; fid < Int(gwt_p.nelements()); ++fid){
444 27 : gwt_p[fid]->get(a_gwt_p);
445 27 : activeFieldIndex_p=fid;
446 27 : if (rmode=="norm") {
447 19 : os << "Normal robustness, robust = " << robust << LogIO::POST;
448 19 : Double sumlocwt = 0.;
449 7707 : for(Int vgrid=0;vgrid<ny;vgrid++) {
450 4252392 : for(Int ugrid=0;ugrid<nx;ugrid++) {
451 4244704 : if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
452 : }
453 : }
454 19 : f2_p[fid] = square(5.0*pow(10.0,Double(-robust))) / (sumlocwt / sumwt[fid]);
455 19 : d2_p[fid] = 1.0;
456 :
457 : }
458 8 : else if (rmode=="abs") {
459 : os << "Absolute robustness, robust = " << robust << ", noise = "
460 2 : << noise.get("Jy").getValue() << "Jy" << LogIO::POST;
461 2 : f2_p[fid] = square(robust);
462 2 : d2_p[fid] = 2.0 * square(noise.get("Jy").getValue());
463 :
464 : }
465 : else {
466 6 : f2_p[fid] = 1.0;
467 6 : d2_p[fid] = 0.0;
468 : }
469 : }
470 18 : }
471 :
472 :
473 2 : VisImagingWeight::VisImagingWeight(ImageInterface<Float>& im) : robust_p(0.0), rmode_p(""), noise_p(Quantity(0.0, "Jy")) {
474 :
475 4 : LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE));
476 :
477 2 : doFilter_p=False;
478 :
479 2 : wgtType_p="uniform";
480 2 : nx_p=im.shape()(0);
481 2 : ny_p=im.shape()(1);
482 2 : DirectionCoordinate dc=im.coordinates().directionCoordinate(0);
483 2 : dc.setWorldAxisUnits(Vector<String>(2, "rad"));
484 2 : Double cellx=dc.increment()(0);
485 2 : Double celly=dc.increment()(1);
486 2 : uscale_p=nx_p*cellx;
487 2 : vscale_p=ny_p*celly;
488 2 : uorigin_p=nx_p/2;
489 2 : vorigin_p=ny_p/2;
490 : //Now to recover from image density and other parameters
491 2 : Int nplanes=1;
492 2 : if(im.shape().nelements()==5)
493 0 : nplanes=im.shape()[4];
494 2 : gwt_p.resize(nplanes, True, False);
495 2 : if(im.shape().nelements()==5){
496 0 : IPosition blc(Vector<Int>(5,0));
497 0 : for (Int fid=0;fid<nplanes;fid++)
498 : {
499 0 : gwt_p[fid]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
500 0 : Array<Float> lala;
501 0 : blc[4]=fid;
502 0 : im.getSlice(lala, blc, IPosition(5, nx_p, ny_p,1,1,1), True);
503 0 : gwt_p[fid]->put( lala.reform(IPosition(2, nx_p, ny_p)));
504 :
505 0 : }
506 0 : }
507 : else{
508 2 : Array<Float> lala;
509 2 : im.get(lala, True);
510 2 : gwt_p[0]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
511 2 : gwt_p[0]->put( lala.reform(IPosition(2, nx_p, ny_p)));
512 :
513 2 : }
514 2 : const TableRecord& rec=im.miscInfo();
515 2 : if(rec.isDefined("d2")){
516 2 : d2_p.resize();
517 2 : rec.get("d2", d2_p);
518 2 : f2_p.resize();
519 2 : rec.get("f2", f2_p);
520 2 : multiFieldMap_p.clear();
521 : Int mapsize;
522 2 : rec.get("multimapsize", mapsize);
523 4 : for(Int k=0; k < mapsize; ++k){
524 2 : String key;
525 : Int val;
526 2 : rec.get("key"+String::toString(k), key);
527 2 : rec.get("val"+String::toString(k), val);
528 : //cerr << "key and id " << key << " "<< val << endl;
529 2 : multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(key, val));
530 2 : }
531 :
532 :
533 : }
534 2 : activeFieldIndex_p=0;
535 2 : gwt_p[0]->get(a_gwt_p);
536 2 : if(rec.isDefined("dofilter")){
537 0 : rec.get("dofilter", doFilter_p);
538 0 : rec.get("rbmaj", rbmaj_p);
539 0 : rec.get("rbmin", rbmin_p);
540 0 : rec.get("cospa", cospa_p);
541 0 : rec.get("sinpa", sinpa_p);
542 : }
543 :
544 :
545 2 : }
546 :
547 :
548 16739 : VisImagingWeight::~VisImagingWeight(){
549 : /*for (uInt fid=0; fid < gwt_p.nelements(); ++fid){
550 : gwt_p[fid].resize();
551 : }*/
552 :
553 16739 : }
554 :
555 2 : Vector<Int> VisImagingWeight::shapeOfdensityGrid(){
556 2 : Vector<Int> retval(3, 0);
557 2 : retval(2)=gwt_p.nelements();
558 2 : if(retval(2) > 0){
559 2 : retval[0]=gwt_p[0]->shape()(0);
560 2 : retval[1]=gwt_p[0]->shape()(1);
561 :
562 : }
563 :
564 2 : return retval;
565 0 : }
566 2 : void VisImagingWeight::toImageInterface(casacore::ImageInterface<casacore::Float>& im){
567 :
568 2 : if( wgtType_p != "uniform")
569 0 : throw(AipsError("cannot save weight density for non-Briggs' weighting schemes"));
570 :
571 2 : IPosition where= IPosition(im.shape().nelements(),0);
572 2 : Int lastAx=where.nelements()-1;
573 4 : for (uInt fid=0;fid<gwt_p.nelements(); ++fid){
574 2 : activeFieldIndex_p=fid;
575 2 : gwt_p[fid]->get(a_gwt_p);
576 2 : where[lastAx]=fid;
577 2 : im.putSlice(a_gwt_p, where);
578 :
579 : }
580 2 : Record rec(im.miscInfo());
581 2 : rec.define("d2", d2_p);
582 2 : rec.define("f2", f2_p);
583 2 : rec.define("numfield", Int(multiFieldMap_p.size()));
584 2 : uInt keycount=0;
585 4 : for (auto iter = multiFieldMap_p.begin( ); iter != multiFieldMap_p.end( ); ++iter, ++keycount){
586 2 : rec.define("key"+String::toString(keycount), iter->first);
587 2 : rec.define("val"+String::toString(keycount), iter->second);
588 : }
589 2 : rec.define("multimapsize",keycount);
590 2 : if(doFilter_p){
591 0 : rec.define("dofilter", doFilter_p);
592 0 : rec.define("rbmaj", rbmaj_p);
593 0 : rec.define("rbmin", rbmin_p);
594 0 : rec.define("cospa", cospa_p);
595 0 : rec.define("sinpa", sinpa_p);
596 : }
597 :
598 :
599 2 : im.setMiscInfo(rec);
600 :
601 2 : }
602 32 : void VisImagingWeight::setFilter(const String& type, const Quantity& bmaj,
603 : const Quantity& bmin, const Quantity& bpa)
604 : {
605 :
606 64 : LogIO os(LogOrigin("VisImagingWeight", "setFilter()", WHERE));
607 :
608 : // Steps : (1) Convert the uvtaper information into "sigma" for the uv-domain Gaussian.
609 : // (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.
610 : // (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.
611 : // (4) Multiply the weightdensity grid (each cell) with the evaluated Gaussian above.
612 :
613 : // There are two ways this input may be supplied : In the image domain or the UV domain.
614 : // 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.
615 : // 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.
616 :
617 32 : if (type=="gaussian") {
618 :
619 : // uvtaper input is supplied in the uv domain as the HWHM of a Gaussian. This is a 'uvdistance' or 'baseline length' interpretation
620 : // The FWHM_uv (in wavelengths) = beam_lambda * 2 to take it from HWHM to FWHM.
621 : // sigma_uv = FWHM_uv / (2 sqrt(2 log2))
622 : // xx = 1 / ( 2 sigma_uv^2) = (log2) / (beam_lambda)^2
623 :
624 32 : Bool lambdafilt=false;
625 :
626 32 : if( bmaj.getUnit().contains("lambda"))
627 8 : lambdafilt=true;
628 32 : if(lambdafilt){
629 : os << "Filtering for Gaussian of shape from read: "
630 16 : << bmaj.get("klambda").getValue() << " by "
631 16 : << bmin.get("klambda").getValue() << " (klambda) at p.a. "
632 24 : << bpa.get("deg").getValue() << " (degrees)" << LogIO::POST;
633 :
634 8 : rbmaj_p=log(2.0)/square(bmaj.get("lambda").getValue());
635 8 : rbmin_p=log(2.0)/square(bmin.get("lambda").getValue());
636 :
637 : }
638 : else{
639 :
640 : // uvtaper input is supplied in the image domain, as the FWHM of a Gaussian. This is an 'convolving' angular resolution interpretation
641 : // This FWHM_lm (in radians) must be converted to a "Sigma" of the Gaussian, and then taken to the UV-domain.
642 : // FWHM_lm = beam_radians
643 : // sigma_lm = FWHM_lm / (2 sqrt(2 log2))
644 : // sigma_uv = 1 / ( 2 pi sigma_lm )
645 : // xx = 1 / ( 2 sigma_uv^2) = ( (pi^2)/(4 log2) ) * (beam_radians)^2
646 :
647 : os << "Filtering for Gaussian of shape: "
648 48 : << bmaj.get("arcsec").getValue() << " by "
649 48 : << bmin.get("arcsec").getValue() << " (arcsec) at p.a. "
650 72 : << bpa.get("deg").getValue() << " (degrees)" << LogIO::POST;
651 :
652 : // Convert to values that we can use
653 24 : Double fact = C::pi*C::pi/(4.0*log(2.0));
654 24 : rbmaj_p = fact*square(bmaj.get("rad").getValue());
655 24 : rbmin_p = fact*square(bmin.get("rad").getValue());
656 : }
657 :
658 32 : Double rbpa = MVAngle(bpa).get("rad").getValue();
659 32 : cospa_p = sin(rbpa);
660 32 : sinpa_p = cos(rbpa);
661 :
662 : os << "Filtering for Gaussian of shape after convention: maj"
663 : << rbmaj_p << " min "
664 : << rbmin_p<< " pa "
665 32 : << rbpa << " " << LogIO::POST;
666 :
667 32 : doFilter_p=true;
668 : }
669 : else {
670 0 : os << "Unknown filtering " << type << LogIO::EXCEPTION;
671 : }
672 :
673 32 : }
674 :
675 :
676 486912 : Bool VisImagingWeight::doFilter() const{
677 :
678 486912 : return doFilter_p;
679 : }
680 :
681 :
682 8640 : void VisImagingWeight::filter(Matrix<Float>& imWeight, const Matrix<Bool>& flag,
683 : const Matrix<Double>& uvw,
684 : const Vector<Double>& frequency, const Matrix<Float>& weight) const{
685 :
686 : // Expecting weight[nchan,nrow], where nchan=1 or nChan (data)
687 : // If nchan=1 (i.e., WEIGHT_SPECTRUM unavailable), then the
688 : // following will be handle the channel degeneracy correctly, using %.
689 :
690 :
691 8640 : Int nRow=imWeight.shape()(1);
692 8640 : Int nChan=imWeight.shape()(0);
693 8640 : Int nChanWt=weight.shape()(0);
694 3041280 : for (Int row=0; row<nRow; row++) {
695 63685440 : for (Int chn=0; chn<nChan; chn++) {
696 60652800 : Double invLambdaC=frequency(chn)/C::c;
697 60652800 : Double u = uvw(0,row);
698 60652800 : Double v = uvw(1,row);
699 60652800 : if(!flag(chn,row) && (weight(chn%nChanWt,row)>0.0) ) {
700 : // 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.
701 : // If the uvtaper Gaussian has positionangle=0, then ru=u, rv=v.
702 : // If the uvtaper Gaussian has positionangle=90, then ru=v, rv= -u
703 57853473 : Double ru = invLambdaC*( cospa_p * u + sinpa_p * v);
704 57853473 : Double rv = invLambdaC*(- sinpa_p * u + cospa_p * v);
705 57853473 : Double filter = exp(-rbmaj_p*square(ru) - rbmin_p*square(rv));
706 57853473 : imWeight(chn,row)*=filter;
707 : }
708 : else {
709 2799327 : imWeight(chn,row)=0.0;
710 : }
711 : }
712 : }
713 :
714 :
715 8640 : }
716 :
717 :
718 15908 : void VisImagingWeight::weightUniform(Matrix<Float>& imWeight, const Matrix<Bool>& flag, const Matrix<Double>& uvw,
719 : const Vector<Double>& frequency,
720 : const Matrix<Float>& weight, const Int msId, const Int fieldId) const{
721 :
722 :
723 : //cout << " WEIG " << nx_p << " " << ny_p << " " << gwt_p[0].shape() << endl;
724 : //cout << "f2 " << f2_p[0] << " d2 " << d2_p[0] << " uscale " << uscale_p << " vscal " << vscale_p << " origs " << uorigin_p << " " << vorigin_p << endl;
725 31816 : String mapid=String::toString(msId)+String("_")+String::toString(fieldId);
726 : //cout << "min max gwt " << min(gwt_p[0]) << " " << max(gwt_p[0]) << " mapid " << mapid <<endl;
727 :
728 15908 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) )
729 0 : throw(AipsError("Imaging weight calculation is requested for a data that was not selected"));
730 :
731 15908 : auto mapvp = multiFieldMap_p.find(mapid);
732 15908 : Int fid = mapvp != multiFieldMap_p.end( ) ? mapvp->second : -1;
733 : //Int ndrop=0;
734 15908 : if(activeFieldIndex_p != fid){
735 44 : a_gwt_p=gwt_p[fid]->get();
736 44 : activeFieldIndex_p=fid;
737 : }
738 15908 : Double sumwt=0.0;
739 15908 : Int nRow=imWeight.shape()(1);
740 15908 : Int nChannel=imWeight.shape()(0);
741 15908 : Int nChanWt=weight.shape()(0);
742 : Float u, v;
743 5304044 : for (Int row=0; row<nRow; row++) {
744 32840712 : for (Int chn=0; chn<nChannel; chn++) {
745 27552576 : if (!flag(chn,row)) {
746 26353578 : Float f=frequency(chn)/C::c;
747 26353578 : u=uvw(0, row)*f;
748 26353578 : v=uvw(1, row)*f;
749 26353578 : Int ucell=Int(uscale_p*u+uorigin_p);
750 26353578 : Int vcell=Int(vscale_p*v+vorigin_p);
751 26353578 : imWeight(chn,row)=weight(chn%nChanWt,row);
752 26353578 : if((ucell>0)&&(ucell<nx_p)&&(vcell>0)&&(vcell<ny_p) &&a_gwt_p(ucell,vcell)>0.0) {
753 26353023 : imWeight(chn,row)/=a_gwt_p(ucell,vcell)*f2_p[fid]+d2_p[fid];
754 26353023 : sumwt+=imWeight(chn,row);
755 :
756 : }
757 : else {
758 555 : imWeight(chn,row)=0.0;
759 : //ndrop++;
760 : }
761 : }
762 : else{
763 1198998 : imWeight(chn,row)=0.0;
764 : }
765 : }
766 : }
767 :
768 15908 : }
769 :
770 : /* unused version?
771 : void VisImagingWeight::weightNatural(Matrix<Float>& imagingWeight, const Matrix<Bool>& flag,
772 : const Matrix<Float>& weight) const{
773 :
774 :
775 :
776 : Int nRow=imagingWeight.shape()(1);
777 : Vector<Float> wgtRow(nRow);
778 :
779 : for (Int row=0; row<nRow; row++) {
780 : wgtRow(row)=max(weight.column(row));
781 : }
782 : weightNatural(imagingWeight, flag, wgtRow);
783 :
784 : }
785 : */
786 463567 : void VisImagingWeight::weightNatural(Matrix<Float>& imagingWeight, const Matrix<Bool>& flag,
787 : const Matrix<Float>& weight) const{
788 :
789 463567 : Double sumwt=0.0;
790 : //cerr << "shape of weight " << weight.shape() << " max " << max (weight.column(0) ) << " max " << min(weight.column(0)) << " " << weight.column(0).shape() << endl;
791 463567 : Int nRow=imagingWeight.shape()(1);
792 463567 : Int nChan=imagingWeight.shape()(0);
793 463567 : Int nChanWt=weight.shape()(0);
794 151441975 : for (Int row=0; row<nRow; row++) {
795 1687459004 : for (Int chn=0; chn<nChan; chn++) {
796 1536480596 : if( !flag(chn,row) ) {
797 1467183490 : imagingWeight(chn,row)=weight(chn%nChanWt,row);
798 1467183490 : sumwt+=imagingWeight(chn,row);
799 : }
800 : else {
801 69297106 : imagingWeight(chn,row)=0.0;
802 : }
803 : }
804 : }
805 :
806 :
807 463567 : }
808 :
809 :
810 1440 : void VisImagingWeight::weightRadial(Matrix<Float>& imagingWeight,
811 : const Matrix<Bool>& flag,
812 : const Matrix<Double>& uvw,
813 : const Vector<Double>& frequency,
814 : const Matrix<Float>& weight) const{
815 :
816 1440 : Double sumwt=0.0;
817 1440 : Int nRow=imagingWeight.shape()(1);
818 1440 : Int nChan=imagingWeight.shape()(0);
819 1440 : Int nChanWt=weight.shape()(0);
820 :
821 506880 : for (Int row=0; row<nRow; row++) {
822 1516320 : for (Int chn=0; chn< nChan; chn++) {
823 1010880 : Float f=frequency(chn)/C::c;
824 1010880 : if( !flag(chn,row) ) {
825 1930716 : imagingWeight(chn,row)=
826 965358 : f*sqrt(square(uvw(0, row))+square(uvw(1, row)))
827 965358 : *weight(chn%nChanWt,row);
828 965358 : sumwt+=imagingWeight(chn,row);
829 : }
830 : else {
831 45522 : imagingWeight(chn,row)=0.0;
832 : }
833 : }
834 : }
835 :
836 :
837 1440 : }
838 :
839 10943 : VisImagingWeight& VisImagingWeight::operator=(const VisImagingWeight& other){
840 10943 : if(this != &other){
841 : // gwt_p=other.gwt_p;
842 :
843 10943 : gwt_p.resize(other.gwt_p.nelements(), true, false);
844 11001 : for (uInt k=0; k < gwt_p.nelements(); ++k){
845 : //gwt_p[k].resize();
846 58 : gwt_p[k]=other.gwt_p[k];
847 : }
848 :
849 :
850 10943 : wgtType_p=other.wgtType_p;
851 10943 : uscale_p=other.uscale_p;
852 10943 : vscale_p=other.vscale_p;
853 10943 : f2_p.resize();
854 10943 : d2_p.resize();
855 10943 : f2_p=other.f2_p;
856 10943 : d2_p=other.d2_p;
857 10943 : uorigin_p=other.uorigin_p;
858 10943 : vorigin_p=other.vorigin_p;
859 10943 : nx_p=other.nx_p;
860 10943 : ny_p=other.ny_p;
861 10943 : doFilter_p=other.doFilter_p;
862 10943 : cospa_p=other.cospa_p;
863 10943 : sinpa_p=other.sinpa_p;
864 10943 : rbmaj_p=other.rbmaj_p;
865 10943 : rbmin_p=other.rbmin_p;
866 10943 : robust_p=other.robust_p;
867 10943 : rmode_p=other.rmode_p;
868 10943 : multiFieldMap_p=other.multiFieldMap_p;
869 : }
870 10943 : return *this;
871 : }
872 :
873 1429209 : String VisImagingWeight::getType() const{
874 :
875 1429209 : return wgtType_p;
876 :
877 : }
878 0 : Bool VisImagingWeight::getWeightDensity (Block<Matrix<Float> >& density){
879 0 : if(wgtType_p != "uniform"){
880 0 : density.resize(0, true, false);
881 0 : return false;
882 : }
883 0 : density.resize(gwt_p.nelements(), true, false);
884 0 : for (uInt k=0; k < gwt_p.nelements(); ++k){
885 0 : density[k].resize(gwt_p[k]->shape());
886 0 : density[k]=gwt_p[k]->get();
887 : }
888 0 : return true;
889 : }
890 0 : void VisImagingWeight::setWeightDensity(const Block<Matrix<Float> >& density){
891 0 : if(wgtType_p=="uniform"){
892 0 : gwt_p.resize(density.nelements(), true, false);
893 0 : for (uInt k=0; k < gwt_p.nelements(); ++k){
894 : //gwt_p[k].resize();
895 0 : gwt_p[k]=new TempLattice<Float>(IPosition(2, density[k].shape()[0], density[k].shape()[1]),0);
896 0 : gwt_p[k]->put(density[k]);
897 : }
898 : //break any old reference
899 0 : f2_p.resize();
900 0 : d2_p.resize();
901 0 : f2_p.resize(gwt_p.nelements());
902 0 : d2_p.resize(gwt_p.nelements());
903 0 : f2_p.set(1.0);
904 0 : d2_p.set(0.0);
905 : //Float f2, d2;
906 0 : for(uInt fid=0; fid < gwt_p.nelements(); ++fid){
907 0 : activeFieldIndex_p=fid;
908 0 : a_gwt_p=gwt_p[fid]->get();
909 0 : if (rmode_p=="norm") {
910 0 : Double sumlocwt = 0.;
911 0 : for(Int vgrid=0;vgrid<a_gwt_p.shape()(1);vgrid++) {
912 0 : for(Int ugrid=0;ugrid<a_gwt_p.shape()(0);ugrid++) {
913 0 : if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
914 : }
915 : }
916 0 : Double sumwt_fid=sum(a_gwt_p);
917 0 : f2_p[fid] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumwt_fid);
918 0 : d2_p[fid] = 1.0;
919 : }
920 0 : else if (rmode_p=="abs") {
921 0 : f2_p[fid] = square(robust_p);
922 0 : d2_p[fid] = 2.0 * square(noise_p.get("Jy").getValue());
923 :
924 : }
925 : else {
926 0 : f2_p[fid] = 1.0;
927 0 : d2_p[fid] = 0.0;
928 : }
929 : }
930 :
931 : }
932 :
933 :
934 0 : }
935 5760 : void VisImagingWeight::cube2Matrix(const Cube<Bool>& fcube, Matrix<Bool>& fMat)
936 : {
937 5760 : fMat.resize(fcube.shape()[1], fcube.shape()[2]);
938 : Bool deleteIt1;
939 : Bool deleteIt2;
940 5760 : const Bool * pcube = fcube.getStorage (deleteIt1);
941 5760 : Bool * pflags = fMat.getStorage (deleteIt2);
942 1958516 : for (uInt row = 0; row < fcube.shape()[2]; row++) {
943 11326992 : for (Int chn = 0; chn < fcube.shape()[1]; chn++) {
944 9374236 : *pflags = *pcube++;
945 18766072 : for (Int pol = 1; pol < fcube.shape()[0]; pol++, pcube++) {
946 9391836 : *pflags = *pcube ? *pcube : *pflags;
947 : }
948 9374236 : pflags++;
949 : }
950 : }
951 5760 : fcube.freeStorage (pcube, deleteIt1);
952 5760 : fMat.putStorage (pflags, deleteIt2);
953 5760 : }
954 :
955 :
956 491746 : void VisImagingWeight::unPolChanWeight(Matrix<Float>& chanRowWt, const Cube<Float>& corrChanRowWt) {
957 :
958 : //cout << "VIW::uPCW" << endl;
959 :
960 491746 : IPosition sh3=corrChanRowWt.shape();
961 491746 : IPosition sh2=sh3.getLast(2);
962 : //cerr << "sh3" << sh3 << " sh2 " << sh2 << endl;
963 491746 : chanRowWt.resize(sh2);
964 : //cout << "assign..." << endl;
965 491746 : chanRowWt=corrChanRowWt(Slice(0,1,1),Slice(),Slice()).reform(sh2);
966 491746 : const Int& nPol=sh3[0];
967 491746 : if (nPol>1) {
968 491746 : chanRowWt += corrChanRowWt(Slice(nPol-1,1,1),Slice(),Slice()).reform(sh2);
969 491746 : chanRowWt /= 2.0f;
970 : }
971 :
972 : // cout << "VIW::uPCW" << endl;
973 :
974 491746 : }
975 :
976 :
977 : }//# NAMESPACE CASA - END
978 :
979 :
|
