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 11 : 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 11 : }
48 :
49 2 : 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 2 : wgtType_p=type;
53 2 : wgtType_p.downcase();
54 2 : if (wgtType_p != "natural" && wgtType_p != "radial"){
55 :
56 0 : throw(AipsError("Programmer error...wrong constructor used"));
57 : }
58 :
59 2 : }
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 0 : 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 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) {
303 :
304 0 : LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE));
305 :
306 :
307 :
308 0 : vi::VisBuffer2* vb=(visIter.getVisBuffer());
309 0 : wgtType_p="uniform";
310 : // Float uscale, vscale;
311 : //Int uorigin, vorigin;
312 0 : Vector<Double> deltas;
313 0 : uscale_p=(nx*cellx.get("rad").getValue());
314 0 : vscale_p=(ny*celly.get("rad").getValue());
315 0 : uorigin_p=nx/2;
316 0 : vorigin_p=ny/2;
317 0 : nx_p=nx;
318 0 : ny_p=ny;
319 : // Simply declare a big matrix
320 : //Matrix<Float> gwt(nx,ny);
321 0 : gwt_p.resize(1);
322 0 : multiFieldMap_p.clear();
323 0 : visIter.originChunks();
324 0 : visIter.origin();
325 0 : String mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]);
326 :
327 0 : multiFieldMap_p.insert( std::pair<casacore::String, casacore::Int>(mapid, 0) );
328 0 : gwt_p[0]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
329 0 : gwt_p[0]->set(0.0);
330 0 : a_gwt_p.resize(nx_p, ny_p);
331 0 : a_gwt_p.set(0.0);
332 :
333 :
334 : // Discover if weightSpectrum non-trivially available
335 0 : Bool doWtSp=visIter.weightSpectrumExists();
336 :
337 0 : Int fields=0;
338 0 : for (visIter.originChunks();visIter.moreChunks();visIter.nextChunk()) {
339 0 : for (visIter.origin();visIter.more();visIter.next()) {
340 0 : if(vb->isNewFieldId()){
341 0 : mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]);
342 0 : if(multiField){
343 0 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) ){
344 0 : fields+=1;
345 0 : gwt_p.resize(fields+1);
346 0 : gwt_p[fields]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
347 0 : gwt_p[fields]->set(0.0);
348 : }
349 : }
350 0 : if( multiFieldMap_p.find(mapid) == multiFieldMap_p.end( ) )
351 0 : multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(mapid, fields));
352 : }
353 : }
354 : }
355 :
356 : Float u, v;
357 0 : Vector<Double> sumwt(fields+1,0.0);
358 0 : f2_p.resize(fields+1);
359 0 : d2_p.resize(fields+1);
360 0 : Int fid=0;
361 0 : for (visIter.originChunks();visIter.moreChunks();visIter.nextChunk()) {
362 0 : for (visIter.origin();visIter.more();visIter.next()) {
363 0 : if(vb->isNewFieldId()){
364 0 : mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]);
365 :
366 0 : if(multiField){
367 0 : if(activeFieldIndex_p >=0)
368 0 : gwt_p[activeFieldIndex_p]->put(a_gwt_p);
369 : }
370 0 : activeFieldIndex_p= multiFieldMap_p.find(mapid) != multiFieldMap_p.end( ) ? multiFieldMap_p[mapid] : -1;
371 0 : if(multiField && activeFieldIndex_p >=0)
372 0 : gwt_p[activeFieldIndex_p]->get(a_gwt_p);
373 : }
374 :
375 0 : auto mapvp = multiFieldMap_p.find(mapid);
376 0 : fid= mapvp != multiFieldMap_p.end( ) ? mapvp->second : -1;
377 0 : Int nRow=vb->nRows();
378 0 : Int nChan=vb->nChannels();
379 :
380 : // Extract weights correctly (use WEIGHT_SPECTRUM, if available)
381 0 : Matrix<Float> wtm;
382 0 : Cube<Float> wtc;
383 0 : if (doWtSp)
384 : // WS available [nCorr,nChan,nRow]
385 0 : wtc.reference(vb->weightSpectrum());
386 : else {
387 0 : Int nCorr=vb->nCorrelations();
388 : // WS UNavailable weight()[nCorr,nRow] --> [nCorr,nChan,nRow]
389 0 : wtc.reference(vb->weight().reform(IPosition(3,nCorr,1,nRow))); // unchan'd weight as single-chan
390 : }
391 0 : unPolChanWeight(wtm,wtc); // Collapse on corr axis
392 :
393 : // Used for mod of chn index below
394 0 : Int nChanWt=wtm.shape()(0);
395 :
396 : //Oww !!! temporary implementation of old vb.flag just to see if things work
397 0 : Matrix<Bool> flag;
398 0 : cube2Matrix(vb->flagCube(), flag);
399 0 : for (Int row=0; row<nRow; row++) {
400 0 : for (Int chn=0; chn<nChan; chn++) {
401 :
402 0 : Float currwt=wtm(chn%nChanWt,row); // the weight for this chan,row
403 :
404 0 : if(!flag(chn,row)) {
405 0 : Float f=vb->getFrequency(row, chn)/C::c;
406 0 : u=vb->uvw()(0,row)*f;
407 0 : v=vb->uvw()(1,row)*f;
408 0 : Int ucell=Int(uscale_p*u+uorigin_p);
409 0 : Int vcell=Int(vscale_p*v+vorigin_p);
410 0 : if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) {
411 0 : for (Int iv=-vBox;iv<=vBox;iv++) {
412 0 : for (Int iu=-uBox;iu<=uBox;iu++) {
413 0 : a_gwt_p(ucell+iu,vcell+iv)+=currwt;
414 0 : sumwt[fid]+=currwt;
415 : }
416 : }
417 : }
418 0 : ucell=Int(-uscale_p*u+uorigin_p);
419 0 : vcell=Int(-vscale_p*v+vorigin_p);
420 0 : if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) {
421 0 : for (Int iv=-vBox;iv<=vBox;iv++) {
422 0 : for (Int iu=-uBox;iu<=uBox;iu++) {
423 0 : a_gwt_p(ucell+iu,vcell+iv)+=currwt;
424 0 : sumwt[fid]+=currwt;
425 : }
426 : }
427 : }
428 : }
429 : }
430 : }
431 0 : }
432 : }
433 :
434 : // make sure last active plane is saved
435 0 : 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 0 : for(fid=0; fid < Int(gwt_p.nelements()); ++fid){
444 0 : gwt_p[fid]->get(a_gwt_p);
445 0 : activeFieldIndex_p=fid;
446 0 : if (rmode=="norm") {
447 0 : os << "Normal robustness, robust = " << robust << LogIO::POST;
448 0 : Double sumlocwt = 0.;
449 0 : for(Int vgrid=0;vgrid<ny;vgrid++) {
450 0 : for(Int ugrid=0;ugrid<nx;ugrid++) {
451 0 : if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid));
452 : }
453 : }
454 0 : f2_p[fid] = square(5.0*pow(10.0,Double(-robust))) / (sumlocwt / sumwt[fid]);
455 0 : d2_p[fid] = 1.0;
456 :
457 : }
458 0 : else if (rmode=="abs") {
459 : os << "Absolute robustness, robust = " << robust << ", noise = "
460 0 : << noise.get("Jy").getValue() << "Jy" << LogIO::POST;
461 0 : f2_p[fid] = square(robust);
462 0 : d2_p[fid] = 2.0 * square(noise.get("Jy").getValue());
463 :
464 : }
465 : else {
466 0 : f2_p[fid] = 1.0;
467 0 : d2_p[fid] = 0.0;
468 : }
469 : }
470 0 : }
471 :
472 :
473 0 : VisImagingWeight::VisImagingWeight(ImageInterface<Float>& im) : robust_p(0.0), rmode_p(""), noise_p(Quantity(0.0, "Jy")) {
474 :
475 0 : LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE));
476 :
477 0 : doFilter_p=False;
478 :
479 0 : wgtType_p="uniform";
480 0 : nx_p=im.shape()(0);
481 0 : ny_p=im.shape()(1);
482 0 : DirectionCoordinate dc=im.coordinates().directionCoordinate(0);
483 0 : dc.setWorldAxisUnits(Vector<String>(2, "rad"));
484 0 : Double cellx=dc.increment()(0);
485 0 : Double celly=dc.increment()(1);
486 0 : uscale_p=nx_p*cellx;
487 0 : vscale_p=ny_p*celly;
488 0 : uorigin_p=nx_p/2;
489 0 : vorigin_p=ny_p/2;
490 : //Now to recover from image density and other parameters
491 0 : Int nplanes=1;
492 0 : if(im.shape().nelements()==5)
493 0 : nplanes=im.shape()[4];
494 0 : gwt_p.resize(nplanes, True, False);
495 0 : 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 0 : Array<Float> lala;
509 0 : im.get(lala, True);
510 0 : gwt_p[0]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0);
511 0 : gwt_p[0]->put( lala.reform(IPosition(2, nx_p, ny_p)));
512 :
513 0 : }
514 0 : const TableRecord& rec=im.miscInfo();
515 0 : if(rec.isDefined("d2")){
516 0 : d2_p.resize();
517 0 : rec.get("d2", d2_p);
518 0 : f2_p.resize();
519 0 : rec.get("f2", f2_p);
520 0 : multiFieldMap_p.clear();
521 : Int mapsize;
522 0 : rec.get("multimapsize", mapsize);
523 0 : for(Int k=0; k < mapsize; ++k){
524 0 : String key;
525 : Int val;
526 0 : rec.get("key"+String::toString(k), key);
527 0 : rec.get("val"+String::toString(k), val);
528 : //cerr << "key and id " << key << " "<< val << endl;
529 0 : multiFieldMap_p.insert(std::pair<casacore::String, casacore::Int>(key, val));
530 0 : }
531 :
532 :
533 : }
534 0 : activeFieldIndex_p=0;
535 0 : gwt_p[0]->get(a_gwt_p);
536 0 : 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 0 : }
546 :
547 :
548 11 : VisImagingWeight::~VisImagingWeight(){
549 : /*for (uInt fid=0; fid < gwt_p.nelements(); ++fid){
550 : gwt_p[fid].resize();
551 : }*/
552 :
553 11 : }
554 :
555 0 : Vector<Int> VisImagingWeight::shapeOfdensityGrid(){
556 0 : Vector<Int> retval(3, 0);
557 0 : retval(2)=gwt_p.nelements();
558 0 : if(retval(2) > 0){
559 0 : retval[0]=gwt_p[0]->shape()(0);
560 0 : retval[1]=gwt_p[0]->shape()(1);
561 :
562 : }
563 :
564 0 : return retval;
565 0 : }
566 0 : void VisImagingWeight::toImageInterface(casacore::ImageInterface<casacore::Float>& im){
567 :
568 0 : if( wgtType_p != "uniform")
569 0 : throw(AipsError("cannot save weight density for non-Briggs' weighting schemes"));
570 :
571 0 : IPosition where= IPosition(im.shape().nelements(),0);
572 0 : Int lastAx=where.nelements()-1;
573 0 : for (uInt fid=0;fid<gwt_p.nelements(); ++fid){
574 0 : activeFieldIndex_p=fid;
575 0 : gwt_p[fid]->get(a_gwt_p);
576 0 : where[lastAx]=fid;
577 0 : im.putSlice(a_gwt_p, where);
578 :
579 : }
580 0 : Record rec(im.miscInfo());
581 0 : rec.define("d2", d2_p);
582 0 : rec.define("f2", f2_p);
583 0 : rec.define("numfield", Int(multiFieldMap_p.size()));
584 0 : uInt keycount=0;
585 0 : for (auto iter = multiFieldMap_p.begin( ); iter != multiFieldMap_p.end( ); ++iter, ++keycount){
586 0 : rec.define("key"+String::toString(keycount), iter->first);
587 0 : rec.define("val"+String::toString(keycount), iter->second);
588 : }
589 0 : rec.define("multimapsize",keycount);
590 0 : 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 0 : im.setMiscInfo(rec);
600 :
601 0 : }
602 0 : void VisImagingWeight::setFilter(const String& type, const Quantity& bmaj,
603 : const Quantity& bmin, const Quantity& bpa)
604 : {
605 :
606 0 : 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 0 : 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 0 : Bool lambdafilt=false;
625 :
626 0 : if( bmaj.getUnit().contains("lambda"))
627 0 : lambdafilt=true;
628 0 : if(lambdafilt){
629 : os << "Filtering for Gaussian of shape from read: "
630 0 : << bmaj.get("klambda").getValue() << " by "
631 0 : << bmin.get("klambda").getValue() << " (klambda) at p.a. "
632 0 : << bpa.get("deg").getValue() << " (degrees)" << LogIO::POST;
633 :
634 0 : rbmaj_p=log(2.0)/square(bmaj.get("lambda").getValue());
635 0 : 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 0 : << bmaj.get("arcsec").getValue() << " by "
649 0 : << bmin.get("arcsec").getValue() << " (arcsec) at p.a. "
650 0 : << bpa.get("deg").getValue() << " (degrees)" << LogIO::POST;
651 :
652 : // Convert to values that we can use
653 0 : Double fact = C::pi*C::pi/(4.0*log(2.0));
654 0 : rbmaj_p = fact*square(bmaj.get("rad").getValue());
655 0 : rbmin_p = fact*square(bmin.get("rad").getValue());
656 : }
657 :
658 0 : Double rbpa = MVAngle(bpa).get("rad").getValue();
659 0 : cospa_p = sin(rbpa);
660 0 : sinpa_p = cos(rbpa);
661 :
662 : os << "Filtering for Gaussian of shape after convention: maj"
663 : << rbmaj_p << " min "
664 : << rbmin_p<< " pa "
665 0 : << rbpa << " " << LogIO::POST;
666 :
667 0 : doFilter_p=true;
668 : }
669 : else {
670 0 : os << "Unknown filtering " << type << LogIO::EXCEPTION;
671 : }
672 :
673 0 : }
674 :
675 :
676 1000 : Bool VisImagingWeight::doFilter() const{
677 :
678 1000 : return doFilter_p;
679 : }
680 :
681 :
682 0 : 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 0 : Int nRow=imWeight.shape()(1);
692 0 : Int nChan=imWeight.shape()(0);
693 0 : Int nChanWt=weight.shape()(0);
694 0 : for (Int row=0; row<nRow; row++) {
695 0 : for (Int chn=0; chn<nChan; chn++) {
696 0 : Double invLambdaC=frequency(chn)/C::c;
697 0 : Double u = uvw(0,row);
698 0 : Double v = uvw(1,row);
699 0 : 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 0 : Double ru = invLambdaC*( cospa_p * u + sinpa_p * v);
704 0 : Double rv = invLambdaC*(- sinpa_p * u + cospa_p * v);
705 0 : Double filter = exp(-rbmaj_p*square(ru) - rbmin_p*square(rv));
706 0 : imWeight(chn,row)*=filter;
707 : }
708 : else {
709 0 : imWeight(chn,row)=0.0;
710 : }
711 : }
712 : }
713 :
714 :
715 0 : }
716 :
717 :
718 0 : 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 0 : 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 0 : 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 0 : auto mapvp = multiFieldMap_p.find(mapid);
732 0 : Int fid = mapvp != multiFieldMap_p.end( ) ? mapvp->second : -1;
733 : //Int ndrop=0;
734 0 : if(activeFieldIndex_p != fid){
735 0 : a_gwt_p=gwt_p[fid]->get();
736 0 : activeFieldIndex_p=fid;
737 : }
738 0 : Double sumwt=0.0;
739 0 : Int nRow=imWeight.shape()(1);
740 0 : Int nChannel=imWeight.shape()(0);
741 0 : Int nChanWt=weight.shape()(0);
742 : Float u, v;
743 0 : for (Int row=0; row<nRow; row++) {
744 0 : for (Int chn=0; chn<nChannel; chn++) {
745 0 : if (!flag(chn,row)) {
746 0 : Float f=frequency(chn)/C::c;
747 0 : u=uvw(0, row)*f;
748 0 : v=uvw(1, row)*f;
749 0 : Int ucell=Int(uscale_p*u+uorigin_p);
750 0 : Int vcell=Int(vscale_p*v+vorigin_p);
751 0 : imWeight(chn,row)=weight(chn%nChanWt,row);
752 0 : if((ucell>0)&&(ucell<nx_p)&&(vcell>0)&&(vcell<ny_p) &&a_gwt_p(ucell,vcell)>0.0) {
753 0 : imWeight(chn,row)/=a_gwt_p(ucell,vcell)*f2_p[fid]+d2_p[fid];
754 0 : sumwt+=imWeight(chn,row);
755 :
756 : }
757 : else {
758 0 : imWeight(chn,row)=0.0;
759 : //ndrop++;
760 : }
761 : }
762 : else{
763 0 : imWeight(chn,row)=0.0;
764 : }
765 : }
766 : }
767 :
768 0 : }
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 800 : void VisImagingWeight::weightNatural(Matrix<Float>& imagingWeight, const Matrix<Bool>& flag,
787 : const Matrix<Float>& weight) const{
788 :
789 800 : 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 800 : Int nRow=imagingWeight.shape()(1);
792 800 : Int nChan=imagingWeight.shape()(0);
793 800 : Int nChanWt=weight.shape()(0);
794 281600 : for (Int row=0; row<nRow; row++) {
795 28360800 : for (Int chn=0; chn<nChan; chn++) {
796 28080000 : if( !flag(chn,row) ) {
797 28080000 : imagingWeight(chn,row)=weight(chn%nChanWt,row);
798 28080000 : sumwt+=imagingWeight(chn,row);
799 : }
800 : else {
801 0 : imagingWeight(chn,row)=0.0;
802 : }
803 : }
804 : }
805 :
806 :
807 800 : }
808 :
809 :
810 0 : 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 0 : Double sumwt=0.0;
817 0 : Int nRow=imagingWeight.shape()(1);
818 0 : Int nChan=imagingWeight.shape()(0);
819 0 : Int nChanWt=weight.shape()(0);
820 :
821 0 : for (Int row=0; row<nRow; row++) {
822 0 : for (Int chn=0; chn< nChan; chn++) {
823 0 : Float f=frequency(chn)/C::c;
824 0 : if( !flag(chn,row) ) {
825 0 : imagingWeight(chn,row)=
826 0 : f*sqrt(square(uvw(0, row))+square(uvw(1, row)))
827 0 : *weight(chn%nChanWt,row);
828 0 : sumwt+=imagingWeight(chn,row);
829 : }
830 : else {
831 0 : imagingWeight(chn,row)=0.0;
832 : }
833 : }
834 : }
835 :
836 :
837 0 : }
838 :
839 6 : VisImagingWeight& VisImagingWeight::operator=(const VisImagingWeight& other){
840 6 : if(this != &other){
841 : // gwt_p=other.gwt_p;
842 :
843 6 : gwt_p.resize(other.gwt_p.nelements(), true, false);
844 6 : for (uInt k=0; k < gwt_p.nelements(); ++k){
845 : //gwt_p[k].resize();
846 0 : gwt_p[k]=other.gwt_p[k];
847 : }
848 :
849 :
850 6 : wgtType_p=other.wgtType_p;
851 6 : uscale_p=other.uscale_p;
852 6 : vscale_p=other.vscale_p;
853 6 : f2_p.resize();
854 6 : d2_p.resize();
855 6 : f2_p=other.f2_p;
856 6 : d2_p=other.d2_p;
857 6 : uorigin_p=other.uorigin_p;
858 6 : vorigin_p=other.vorigin_p;
859 6 : nx_p=other.nx_p;
860 6 : ny_p=other.ny_p;
861 6 : doFilter_p=other.doFilter_p;
862 6 : cospa_p=other.cospa_p;
863 6 : sinpa_p=other.sinpa_p;
864 6 : rbmaj_p=other.rbmaj_p;
865 6 : rbmin_p=other.rbmin_p;
866 6 : robust_p=other.robust_p;
867 6 : rmode_p=other.rmode_p;
868 6 : multiFieldMap_p=other.multiFieldMap_p;
869 : }
870 6 : return *this;
871 : }
872 :
873 2600 : String VisImagingWeight::getType() const{
874 :
875 2600 : 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 0 : void VisImagingWeight::cube2Matrix(const Cube<Bool>& fcube, Matrix<Bool>& fMat)
936 : {
937 0 : fMat.resize(fcube.shape()[1], fcube.shape()[2]);
938 : Bool deleteIt1;
939 : Bool deleteIt2;
940 0 : const Bool * pcube = fcube.getStorage (deleteIt1);
941 0 : Bool * pflags = fMat.getStorage (deleteIt2);
942 0 : for (uInt row = 0; row < fcube.shape()[2]; row++) {
943 0 : for (Int chn = 0; chn < fcube.shape()[1]; chn++) {
944 0 : *pflags = *pcube++;
945 0 : for (Int pol = 1; pol < fcube.shape()[0]; pol++, pcube++) {
946 0 : *pflags = *pcube ? *pcube : *pflags;
947 : }
948 0 : pflags++;
949 : }
950 : }
951 0 : fcube.freeStorage (pcube, deleteIt1);
952 0 : fMat.putStorage (pflags, deleteIt2);
953 0 : }
954 :
955 :
956 600 : void VisImagingWeight::unPolChanWeight(Matrix<Float>& chanRowWt, const Cube<Float>& corrChanRowWt) {
957 :
958 : //cout << "VIW::uPCW" << endl;
959 :
960 600 : IPosition sh3=corrChanRowWt.shape();
961 600 : IPosition sh2=sh3.getLast(2);
962 : //cerr << "sh3" << sh3 << " sh2 " << sh2 << endl;
963 600 : chanRowWt.resize(sh2);
964 : //cout << "assign..." << endl;
965 600 : chanRowWt=corrChanRowWt(Slice(0,1,1),Slice(),Slice()).reform(sh2);
966 600 : const Int& nPol=sh3[0];
967 600 : if (nPol>1) {
968 600 : chanRowWt += corrChanRowWt(Slice(nPol-1,1,1),Slice(),Slice()).reform(sh2);
969 600 : chanRowWt /= 2.0f;
970 : }
971 :
972 : // cout << "VIW::uPCW" << endl;
973 :
974 600 : }
975 :
976 :
977 : }//# NAMESPACE CASA - END
978 :
979 :
|
