Commits
Ville Suoranta authored e653b02bfe1 Merge
1 1 | //# VisImagingWeight.cc: imaging weight calculation for a give buffer |
2 - | //# Copyright (C) 2009 |
2 + | //# Copyright (C) 2009-2018 |
3 3 | //# Associated Universities, Inc. Washington DC, USA. |
4 4 | //# |
5 5 | //# This library is free software; you can redistribute it and/or modify it |
6 - | //# under the terms of the GNU Library General Public License as published by |
6 + | //# under the terms of the GNU General Public License as published by |
7 7 | //# the Free Software Foundation; either version 2 of the License, or (at your |
8 8 | //# option) any later version. |
9 9 | //# |
10 10 | //# This library is distributed in the hope that it will be useful, but WITHOUT |
11 11 | //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
12 - | //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public |
12 + | //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
13 13 | //# License for more details. |
14 14 | //# |
15 - | //# You should have received a copy of the GNU Library General Public License |
15 + | //# You should have received a copy of the GNU General Public License |
16 16 | //# along with this library; if not, write to the Free Software Foundation, |
17 17 | //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. |
18 18 | //# |
19 19 | //# Correspondence concerning AIPS++ should be adressed as follows: |
20 20 | //# Internet email: aips2-request@nrao.edu. |
21 21 | //# Postal address: AIPS++ Project Office |
22 22 | //# National Radio Astronomy Observatory |
23 23 | //# 520 Edgemont Road |
24 24 | //# Charlottesville, VA 22903-2475 USA |
25 25 | //# |
26 26 | //# |
27 27 | //# $Id$ |
28 28 | |
29 29 | |
30 30 | |
31 31 | |
32 32 | |
33 33 | |
34 34 | |
35 35 | |
36 - | |
36 + | |
37 37 | |
38 38 | |
39 39 | |
40 40 | |
41 - | |
42 41 | using namespace casacore; |
43 42 | namespace casa { //# NAMESPACE CASA - BEGIN |
44 43 | |
45 - | VisImagingWeight::VisImagingWeight() : multiFieldMap_p(-1), wgtType_p("none"), doFilter_p(false), robust_p(0.0), rmode_p("norm"), noise_p(Quantity(0.0, "Jy")) { |
44 + | VisImagingWeight::VisImagingWeight() : multiFieldMap_p(-1), wgtType_p("none"), doFilter_p(false), robust_p(0.0), rmode_p("norm"), noise_p(Quantity(0.0, "Jy")), activeFieldIndex_p(-1) { |
46 45 | |
47 46 | } |
48 47 | |
49 - | VisImagingWeight::VisImagingWeight(const String& type) : multiFieldMap_p(-1),doFilter_p(false), robust_p(0.0), rmode_p("norm"), noise_p(Quantity(0.0, "Jy")) { |
48 + | VisImagingWeight::VisImagingWeight(const String& type) : multiFieldMap_p(-1),doFilter_p(false), robust_p(0.0), rmode_p("norm"), noise_p(Quantity(0.0, "Jy")), activeFieldIndex_p(-1) { |
50 49 | |
51 50 | wgtType_p=type; |
52 51 | wgtType_p.downcase(); |
53 52 | if (wgtType_p != "natural" && wgtType_p != "radial"){ |
54 53 | |
55 54 | throw(AipsError("Programmer error...wrong constructor used")); |
56 55 | } |
57 56 | |
58 57 | } |
59 58 | |
60 59 | |
61 60 | VisImagingWeight::VisImagingWeight(ROVisibilityIterator& vi, const String& rmode, const Quantity& noise, |
62 61 | const Double robust, const Int nx, const Int ny, |
63 62 | const Quantity& cellx, const Quantity& celly, |
64 - | const Int uBox, const Int vBox, const Bool multiField) : multiFieldMap_p(-1), doFilter_p(false), robust_p(robust), rmode_p(rmode), noise_p(noise) { |
63 + | const Int uBox, const Int vBox, const Bool multiField) : multiFieldMap_p(-1), doFilter_p(false), robust_p(robust), rmode_p(rmode), noise_p(noise), activeFieldIndex_p(-1) { |
65 64 | |
66 65 | LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE)); |
67 66 | |
68 67 | VisBufferAutoPtr vb (vi); |
69 68 | wgtType_p="uniform"; |
70 69 | // Float uscale, vscale; |
71 70 | //Int uorigin, vorigin; |
72 71 | Vector<Double> deltas; |
73 72 | uscale_p=(nx*cellx.get("rad").getValue()); |
74 73 | vscale_p=(ny*celly.get("rad").getValue()); |
81 80 | gwt_p.resize(1); |
82 81 | multiFieldMap_p.clear(); |
83 82 | vi.originChunks(); |
84 83 | vi.origin(); |
85 84 | |
86 85 | // Detect usable WEIGHT_SPECTRUM |
87 86 | Bool doWtSp=(vb->weightSpectrum().nelements()>0); |
88 87 | |
89 88 | String mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()); |
90 89 | multiFieldMap_p.define(mapid, 0); |
91 - | gwt_p[0].resize(nx, ny); |
92 - | gwt_p[0].set(0.0); |
90 + | gwt_p[0]=new TempLattice<Float>(IPosition(2, nx,ny), 0); |
91 + | gwt_p[0]->set(0.0); |
92 + | a_gwt_p.resize(nx_p, ny_p); |
93 + | a_gwt_p.set(0.0); |
93 94 | |
94 95 | Int fields=0; |
95 96 | for (vi.originChunks();vi.moreChunks();vi.nextChunk()) { |
96 97 | for (vi.origin();vi.more();vi++) { |
98 + | |
97 99 | if(vb->newFieldId()){ |
98 100 | mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()); |
99 101 | if(multiField){ |
100 102 | if(!multiFieldMap_p.isDefined(mapid)){ |
101 103 | fields+=1; |
102 104 | gwt_p.resize(fields+1); |
103 - | gwt_p[fields].resize(nx,ny); |
104 - | gwt_p[fields].set(0.0); |
105 + | gwt_p[fields]=new TempLattice<Float>(IPosition(2, nx_p,ny_p),0); |
106 + | gwt_p[fields]->set(0.0); |
105 107 | } |
106 108 | } |
107 109 | if(!multiFieldMap_p.isDefined(mapid)) |
108 110 | multiFieldMap_p.define(mapid, fields); |
111 + | |
109 112 | } |
110 113 | } |
111 114 | } |
112 115 | |
113 116 | Float u, v; |
114 117 | Vector<Double> sumwt(fields+1,0.0); |
115 118 | f2_p.resize(fields+1); |
116 119 | d2_p.resize(fields+1); |
117 120 | Int fid=0; |
118 121 | for (vi.originChunks();vi.moreChunks();vi.nextChunk()) { |
119 122 | for (vi.origin();vi.more();vi++) { |
120 - | if(vb->newFieldId()) |
123 + | if(vb->newFieldId()){ |
121 124 | mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()); |
125 + | if(multiField){ |
126 + | if(activeFieldIndex_p >=0) |
127 + | gwt_p[activeFieldIndex_p]->put(a_gwt_p); |
128 + | } |
129 + | activeFieldIndex_p=multiFieldMap_p(mapid); |
130 + | if(multiField && activeFieldIndex_p >=0) |
131 + | gwt_p[activeFieldIndex_p]->get(a_gwt_p); |
132 + | |
133 + | } |
122 134 | fid=multiFieldMap_p(mapid); |
123 135 | Int nRow=vb->nRow(); |
124 136 | Int nChan=vb->nChannel(); |
125 137 | |
126 138 | // Extract weights correctly (use WEIGHT_SPECTRUM, if available) |
127 139 | Matrix<Float> wtm; |
128 140 | //////////////// |
129 141 | doWtSp=(vb->weightSpectrum().nelements()>0); |
130 142 | ////////////// |
131 143 | if (doWtSp) |
144 156 | Float currwt=wtm(chn%nChanWt,row); // the weight for this chan,row |
145 157 | Float f=vb->frequency()(chn)/C::c; |
146 158 | u=vb->uvw()(row)(0)*f; |
147 159 | v=vb->uvw()(row)(1)*f; |
148 160 | Int ucell=Int(uscale_p*u+uorigin_p); |
149 161 | Int vcell=Int(vscale_p*v+vorigin_p); |
150 162 | |
151 163 | if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) { |
152 164 | for (Int iv=-vBox;iv<=vBox;iv++) { |
153 165 | for (Int iu=-uBox;iu<=uBox;iu++) { |
154 - | gwt_p[fid](ucell+iu,vcell+iv)+=currwt; |
166 + | a_gwt_p(ucell+iu,vcell+iv)+=currwt; |
155 167 | sumwt[fid]+=currwt; |
156 168 | } |
157 169 | } |
158 170 | } |
159 171 | ucell=Int(-uscale_p*u+uorigin_p); |
160 172 | vcell=Int(-vscale_p*v+vorigin_p); |
161 173 | if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) { |
162 174 | for (Int iv=-vBox;iv<=vBox;iv++) { |
163 175 | for (Int iu=-uBox;iu<=uBox;iu++) { |
164 - | gwt_p[fid](ucell+iu,vcell+iv)+=currwt; |
176 + | a_gwt_p(ucell+iu,vcell+iv)+=currwt; |
165 177 | sumwt[fid]+=currwt; |
166 178 | } |
167 179 | } |
168 180 | } |
169 181 | } |
170 182 | } |
171 183 | } |
172 184 | } |
173 185 | } |
174 - | |
186 + | // make sure last active plane is saved |
187 + | gwt_p[activeFieldIndex_p]->put(a_gwt_p); |
175 188 | // We use the approximation that all statistical weights are equal to |
176 189 | // calculate the average summed weights (over visibilities, not bins!) |
177 190 | // This is simply to try an ensure that the normalization of the robustness |
178 191 | // parameter is similar to that of the ungridded case, but it doesn't have |
179 192 | // to be exact, since any given case will require some experimentation. |
180 193 | |
181 194 | //Float f2, d2; |
182 195 | for(fid=0; fid < Int(gwt_p.nelements()); ++fid){ |
196 + | gwt_p[fid]->get(a_gwt_p); |
197 + | activeFieldIndex_p=fid; |
183 198 | if (rmode=="norm") { |
184 199 | os << "Normal robustness, robust = " << robust << LogIO::POST; |
185 200 | Double sumlocwt = 0.; |
186 201 | for(Int vgrid=0;vgrid<ny;vgrid++) { |
187 202 | for(Int ugrid=0;ugrid<nx;ugrid++) { |
188 - | if(gwt_p[fid](ugrid, vgrid)>0.0) sumlocwt+=square(gwt_p[fid](ugrid,vgrid)); |
203 + | if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid)); |
189 204 | } |
190 205 | } |
191 206 | f2_p[fid] = square(5.0*pow(10.0,Double(-robust))) / (sumlocwt / sumwt[fid]); |
192 207 | d2_p[fid] = 1.0; |
193 208 | |
194 209 | } |
195 210 | else if (rmode=="abs") { |
196 211 | os << "Absolute robustness, robust = " << robust << ", noise = " |
197 212 | << noise.get("Jy").getValue() << "Jy" << LogIO::POST; |
198 213 | f2_p[fid] = square(robust); |
199 214 | d2_p[fid] = 2.0 * square(noise.get("Jy").getValue()); |
200 215 | |
201 216 | } |
202 217 | else { |
203 218 | f2_p[fid] = 1.0; |
204 219 | d2_p[fid] = 0.0; |
205 220 | } |
206 221 | } |
207 222 | |
208 223 | } |
209 224 | |
210 - | VisImagingWeight::VisImagingWeight(ROVisibilityIterator& vi, Block<Matrix<Float> >& grids, const String& rmode, const Quantity& noise, |
225 + | VisImagingWeight::VisImagingWeight(ROVisibilityIterator& vi, Block<CountedPtr<TempLattice<Float> > >& grids, const String& rmode, const Quantity& noise, |
211 226 | const Double robust, const Quantity& cellx, const Quantity& celly, |
212 227 | const Bool multiField) : multiFieldMap_p(-1), doFilter_p(false), robust_p(robust), rmode_p(rmode), noise_p(noise) { |
213 228 | |
214 229 | LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE)); |
215 230 | |
216 231 | VisBufferAutoPtr vb (vi); |
217 232 | wgtType_p="uniform"; |
218 233 | gwt_p.resize(grids.nelements(), true, false); |
219 234 | for (uInt k =0; k < gwt_p.nelements(); ++k) |
220 - | gwt_p[k].assign(grids[k]); |
221 - | nx_p=gwt_p[0].shape()[0]; |
222 - | ny_p=gwt_p[0].shape()[1]; |
235 + | gwt_p[k]=grids[k]; |
236 + | nx_p=gwt_p[0]->shape()[0]; |
237 + | ny_p=gwt_p[0]->shape()[1]; |
223 238 | uscale_p=(nx_p*cellx.get("rad").getValue()); |
224 239 | vscale_p=(ny_p*celly.get("rad").getValue()); |
225 240 | uorigin_p=nx_p/2; |
226 241 | vorigin_p=ny_p/2; |
227 242 | |
228 243 | multiFieldMap_p.clear(); |
229 244 | Int fields=0; |
230 245 | String mapid; |
231 246 | for (vi.originChunks();vi.moreChunks();vi.nextChunk()) { |
232 247 | for (vi.origin();vi.more();vi++) { |
239 254 | } |
240 255 | } |
241 256 | if(!multiFieldMap_p.isDefined(mapid)) |
242 257 | multiFieldMap_p.define(mapid, fields); |
243 258 | } |
244 259 | } |
245 260 | } |
246 261 | f2_p.resize(gwt_p.nelements()); |
247 262 | d2_p.resize(gwt_p.nelements()); |
248 263 | for(uInt fid=0; fid < gwt_p.nelements(); ++fid){ |
264 + | activeFieldIndex_p=fid; |
265 + | gwt_p[fid]->get(a_gwt_p); |
249 266 | if (rmode_p=="norm") { |
250 267 | os << "Normal robustness, robust = " << robust_p << LogIO::POST; |
251 268 | Double sumlocwt = 0.; |
252 269 | for(Int vgrid=0;vgrid<ny_p;vgrid++) { |
253 270 | for(Int ugrid=0;ugrid<nx_p;ugrid++) { |
254 - | if(gwt_p[fid](ugrid, vgrid)>0.0) sumlocwt+=square(gwt_p[fid](ugrid,vgrid)); |
271 + | if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid)); |
255 272 | } |
256 273 | } |
257 - | Double sumwt_fid=sum(gwt_p[fid]); |
274 + | Double sumwt_fid=sum(a_gwt_p); |
258 275 | f2_p[fid] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumwt_fid); |
259 276 | d2_p[fid] = 1.0; |
260 277 | |
261 278 | } |
262 279 | else if (rmode=="abs") { |
263 280 | os << "Absolute robustness, robust = " << robust_p << ", noise = " |
264 281 | << noise_p.get("Jy").getValue() << "Jy" << LogIO::POST; |
265 282 | f2_p[fid] = square(robust_p); |
266 283 | d2_p[fid] = 2.0 * square(noise_p.get("Jy").getValue()); |
267 284 | |
270 287 | f2_p[fid] = 1.0; |
271 288 | d2_p[fid] = 0.0; |
272 289 | } |
273 290 | } |
274 291 | } |
275 292 | |
276 293 | VisImagingWeight::VisImagingWeight(vi::VisibilityIterator2& visIter, |
277 294 | const String& rmode, const Quantity& noise, |
278 295 | const Double robust, const Int nx, const Int ny, |
279 296 | const Quantity& cellx, const Quantity& celly, |
280 - | const Int uBox, const Int vBox, const Bool multiField) : multiFieldMap_p(-1), doFilter_p(false), robust_p(robust), rmode_p(rmode), noise_p(noise) { |
297 + | const Int uBox, const Int vBox, const Bool multiField) : multiFieldMap_p(-1), doFilter_p(false), robust_p(robust), rmode_p(rmode), noise_p(noise), activeFieldIndex_p(-1) { |
281 298 | |
282 299 | LogIO os(LogOrigin("VisSetUtil", "VisImagingWeight()", WHERE)); |
283 300 | |
284 301 | |
285 302 | |
286 303 | vi::VisBuffer2* vb=(visIter.getVisBuffer()); |
287 304 | wgtType_p="uniform"; |
288 305 | // Float uscale, vscale; |
289 306 | //Int uorigin, vorigin; |
290 307 | Vector<Double> deltas; |
295 312 | nx_p=nx; |
296 313 | ny_p=ny; |
297 314 | // Simply declare a big matrix |
298 315 | //Matrix<Float> gwt(nx,ny); |
299 316 | gwt_p.resize(1); |
300 317 | multiFieldMap_p.clear(); |
301 318 | visIter.originChunks(); |
302 319 | visIter.origin(); |
303 320 | String mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]); |
304 321 | multiFieldMap_p.define(mapid, 0); |
305 - | gwt_p[0].resize(nx, ny); |
306 - | gwt_p[0].set(0.0); |
322 + | gwt_p[0]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0); |
323 + | gwt_p[0]->set(0.0); |
324 + | a_gwt_p.resize(nx_p, ny_p); |
325 + | a_gwt_p.set(0.0); |
307 326 | |
308 327 | // Discover if weightSpectrum non-trivially available |
309 328 | Bool doWtSp=visIter.weightSpectrumExists(); |
310 329 | |
311 330 | Int fields=0; |
312 331 | for (visIter.originChunks();visIter.moreChunks();visIter.nextChunk()) { |
313 332 | for (visIter.origin();visIter.more();visIter.next()) { |
314 333 | if(vb->isNewFieldId()){ |
315 334 | mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]); |
316 335 | if(multiField){ |
317 336 | if(!multiFieldMap_p.isDefined(mapid)){ |
318 337 | fields+=1; |
319 338 | gwt_p.resize(fields+1); |
320 - | gwt_p[fields].resize(nx,ny); |
321 - | gwt_p[fields].set(0.0); |
339 + | gwt_p[fields]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0); |
340 + | gwt_p[fields]->set(0.0); |
322 341 | } |
323 342 | } |
324 343 | if(!multiFieldMap_p.isDefined(mapid)) |
325 344 | multiFieldMap_p.define(mapid, fields); |
326 345 | } |
327 346 | } |
328 347 | } |
329 348 | |
330 349 | Float u, v; |
331 350 | Vector<Double> sumwt(fields+1,0.0); |
332 351 | f2_p.resize(fields+1); |
333 352 | d2_p.resize(fields+1); |
334 353 | Int fid=0; |
335 354 | for (visIter.originChunks();visIter.moreChunks();visIter.nextChunk()) { |
336 355 | for (visIter.origin();visIter.more();visIter.next()) { |
337 - | if(vb->isNewFieldId()) |
356 + | if(vb->isNewFieldId()){ |
338 357 | mapid=String::toString(vb->msId())+String("_")+String::toString(vb->fieldId()[0]); |
358 + | if(multiField){ |
359 + | if(activeFieldIndex_p >=0) |
360 + | gwt_p[activeFieldIndex_p]->put(a_gwt_p); |
361 + | } |
362 + | activeFieldIndex_p=multiFieldMap_p(mapid); |
363 + | if(multiField && activeFieldIndex_p >=0) |
364 + | gwt_p[activeFieldIndex_p]->get(a_gwt_p); |
365 + | } |
339 366 | fid=multiFieldMap_p(mapid); |
340 367 | Int nRow=vb->nRows(); |
341 368 | Int nChan=vb->nChannels(); |
342 369 | |
343 370 | // Extract weights correctly (use WEIGHT_SPECTRUM, if available) |
344 371 | Matrix<Float> wtm; |
345 372 | Cube<Float> wtc; |
346 373 | if (doWtSp) |
347 374 | // WS available [nCorr,nChan,nRow] |
348 375 | wtc.reference(vb->weightSpectrum()); |
366 393 | |
367 394 | if(!flag(chn,row)) { |
368 395 | Float f=vb->getFrequency(row, chn)/C::c; |
369 396 | u=vb->uvw()(0,row)*f; |
370 397 | v=vb->uvw()(1,row)*f; |
371 398 | Int ucell=Int(uscale_p*u+uorigin_p); |
372 399 | Int vcell=Int(vscale_p*v+vorigin_p); |
373 400 | if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) { |
374 401 | for (Int iv=-vBox;iv<=vBox;iv++) { |
375 402 | for (Int iu=-uBox;iu<=uBox;iu++) { |
376 - | gwt_p[fid](ucell+iu,vcell+iv)+=currwt; |
403 + | a_gwt_p(ucell+iu,vcell+iv)+=currwt; |
377 404 | sumwt[fid]+=currwt; |
378 405 | } |
379 406 | } |
380 407 | } |
381 408 | ucell=Int(-uscale_p*u+uorigin_p); |
382 409 | vcell=Int(-vscale_p*v+vorigin_p); |
383 410 | if(((ucell-uBox)>0)&&((ucell+uBox)<nx)&&((vcell-vBox)>0)&&((vcell+vBox)<ny)) { |
384 411 | for (Int iv=-vBox;iv<=vBox;iv++) { |
385 412 | for (Int iu=-uBox;iu<=uBox;iu++) { |
386 - | gwt_p[fid](ucell+iu,vcell+iv)+=currwt; |
413 + | a_gwt_p(ucell+iu,vcell+iv)+=currwt; |
387 414 | sumwt[fid]+=currwt; |
388 415 | } |
389 416 | } |
390 417 | } |
391 418 | } |
392 419 | } |
393 420 | } |
394 421 | } |
395 422 | } |
396 423 | |
424 + | // make sure last active plane is saved |
425 + | gwt_p[activeFieldIndex_p]->put(a_gwt_p); |
397 426 | // We use the approximation that all statistical weights are equal to |
398 427 | // calculate the average summed weights (over visibilities, not bins!) |
399 428 | // This is simply to try an ensure that the normalization of the robustness |
400 429 | // parameter is similar to that of the ungridded case, but it doesn't have |
401 430 | // to be exact, since any given case will require some experimentation. |
402 431 | |
403 432 | //Float f2, d2; |
404 433 | for(fid=0; fid < Int(gwt_p.nelements()); ++fid){ |
434 + | gwt_p[fid]->get(a_gwt_p); |
435 + | activeFieldIndex_p=fid; |
405 436 | if (rmode=="norm") { |
406 437 | os << "Normal robustness, robust = " << robust << LogIO::POST; |
407 438 | Double sumlocwt = 0.; |
408 439 | for(Int vgrid=0;vgrid<ny;vgrid++) { |
409 440 | for(Int ugrid=0;ugrid<nx;ugrid++) { |
410 - | if(gwt_p[fid](ugrid, vgrid)>0.0) sumlocwt+=square(gwt_p[fid](ugrid,vgrid)); |
441 + | if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid)); |
411 442 | } |
412 443 | } |
413 444 | f2_p[fid] = square(5.0*pow(10.0,Double(-robust))) / (sumlocwt / sumwt[fid]); |
414 445 | d2_p[fid] = 1.0; |
415 446 | |
416 447 | } |
417 448 | else if (rmode=="abs") { |
418 449 | os << "Absolute robustness, robust = " << robust << ", noise = " |
419 450 | << noise.get("Jy").getValue() << "Jy" << LogIO::POST; |
420 451 | f2_p[fid] = square(robust); |
448 479 | vorigin_p=ny_p/2; |
449 480 | //Now to recover from image density and other parameters |
450 481 | Int nplanes=1; |
451 482 | if(im.shape().nelements()==5) |
452 483 | nplanes=im.shape()[4]; |
453 484 | gwt_p.resize(nplanes, True, False); |
454 485 | if(im.shape().nelements()==5){ |
455 486 | IPosition blc(Vector<Int>(5,0)); |
456 487 | for (Int fid=0;fid<nplanes;fid++) |
457 488 | { |
458 - | gwt_p[fid].resize(); |
489 + | gwt_p[fid]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0); |
459 490 | Array<Float> lala; |
460 491 | blc[4]=fid; |
461 492 | im.getSlice(lala, blc, IPosition(5, nx_p, ny_p,1,1,1), True); |
462 - | gwt_p[fid].reference( lala.reform(IPosition(2, nx_p, ny_p))); |
493 + | gwt_p[fid]->put( lala.reform(IPosition(2, nx_p, ny_p))); |
494 + | |
463 495 | } |
464 496 | } |
465 497 | else{ |
466 498 | Array<Float> lala; |
467 499 | im.get(lala, True); |
468 - | gwt_p[0].reference( lala.reform(IPosition(2, nx_p, ny_p))); |
500 + | gwt_p[0]=new TempLattice<Float>(IPosition(2,nx_p, ny_p), 0); |
501 + | gwt_p[0]->put( lala.reform(IPosition(2, nx_p, ny_p))); |
502 + | |
469 503 | } |
470 504 | const TableRecord& rec=im.miscInfo(); |
471 505 | if(rec.isDefined("d2")){ |
472 506 | d2_p.resize(); |
473 507 | rec.get("d2", d2_p); |
474 508 | f2_p.resize(); |
475 509 | rec.get("f2", f2_p); |
476 510 | multiFieldMap_p.clear(); |
477 511 | for(Int k=0; k < nplanes; ++k){ |
478 512 | String key; |
479 513 | Int val; |
480 514 | rec.get("key"+String::toString(k), key); |
481 515 | rec.get("val"+String::toString(k), val); |
482 516 | multiFieldMap_p.define(key, val); |
483 517 | } |
484 518 | |
485 519 | |
486 520 | } |
521 + | activeFieldIndex_p=0; |
522 + | gwt_p[0]->get(a_gwt_p); |
523 + | |
487 524 | |
488 525 | } |
489 526 | |
490 527 | |
491 528 | VisImagingWeight::~VisImagingWeight(){ |
492 - | for (uInt fid=0; fid < gwt_p.nelements(); ++fid){ |
529 + | /*for (uInt fid=0; fid < gwt_p.nelements(); ++fid){ |
493 530 | gwt_p[fid].resize(); |
494 - | } |
531 + | }*/ |
532 + | |
495 533 | } |
496 534 | |
497 535 | Vector<Int> VisImagingWeight::shapeOfdensityGrid(){ |
498 536 | Vector<Int> retval(3, 0); |
499 537 | retval(2)=gwt_p.nelements(); |
500 538 | if(retval(2) > 0){ |
501 - | retval[0]=gwt_p[0].shape()(0); |
502 - | retval[1]=gwt_p[0].shape()(1); |
539 + | retval[0]=gwt_p[0]->shape()(0); |
540 + | retval[1]=gwt_p[0]->shape()(1); |
541 + | |
503 542 | } |
504 543 | |
505 544 | return retval; |
506 545 | } |
507 546 | void VisImagingWeight::toImageInterface(casacore::ImageInterface<casacore::Float>& im){ |
508 547 | |
509 548 | if( wgtType_p != "uniform") |
510 549 | throw(AipsError("cannot save weight density for non-Briggs' weighting schemes")); |
511 550 | |
512 551 | IPosition where= IPosition(im.shape().nelements(),0); |
513 552 | Int lastAx=where.nelements()-1; |
514 553 | for (uInt fid=0;fid<gwt_p.nelements(); ++fid){ |
554 + | activeFieldIndex_p=fid; |
555 + | gwt_p[fid]->get(a_gwt_p); |
515 556 | where[lastAx]=fid; |
516 - | im.putSlice(gwt_p[fid], where); |
557 + | im.putSlice(a_gwt_p, where); |
558 + | |
517 559 | } |
518 560 | Record rec; |
519 561 | rec.define("d2", d2_p); |
520 562 | rec.define("f2", f2_p); |
521 563 | rec.define("numfield", Int(multiFieldMap_p.ndefined())); |
522 564 | for(uInt k=0; k < multiFieldMap_p.ndefined(); ++k){ |
523 565 | rec.define("key"+String::toString(k), multiFieldMap_p.getKey(k)); |
524 566 | rec.define("val"+String::toString(k), multiFieldMap_p.getVal(k)); |
525 567 | } |
526 568 | im.setMiscInfo(rec); |
619 661 | |
620 662 | //cout << " WEIG " << nx_p << " " << ny_p << " " << gwt_p[0].shape() << endl; |
621 663 | //cout << "f2 " << f2_p[0] << " d2 " << d2_p[0] << " uscale " << uscale_p << " vscal " << vscale_p << " origs " << uorigin_p << " " << vorigin_p << endl; |
622 664 | String mapid=String::toString(msId)+String("_")+String::toString(fieldId); |
623 665 | //cout << "min max gwt " << min(gwt_p[0]) << " " << max(gwt_p[0]) << " mapid " << mapid <<endl; |
624 666 | if(!multiFieldMap_p.isDefined(mapid)) |
625 667 | throw(AipsError("Imaging weight calculation is requested for a data that was not selected")); |
626 668 | |
627 669 | Int fid=multiFieldMap_p(mapid); |
628 670 | //Int ndrop=0; |
629 - | |
671 + | if(activeFieldIndex_p != fid){ |
672 + | a_gwt_p=gwt_p[fid]->get(); |
673 + | activeFieldIndex_p=fid; |
674 + | } |
630 675 | Double sumwt=0.0; |
631 676 | Int nRow=imWeight.shape()(1); |
632 677 | Int nChannel=imWeight.shape()(0); |
633 678 | Int nChanWt=weight.shape()(0); |
634 679 | Float u, v; |
635 680 | for (Int row=0; row<nRow; row++) { |
636 681 | for (Int chn=0; chn<nChannel; chn++) { |
637 682 | if (!flag(chn,row)) { |
638 683 | Float f=frequency(chn)/C::c; |
639 684 | u=uvw(0, row)*f; |
640 685 | v=uvw(1, row)*f; |
641 686 | Int ucell=Int(uscale_p*u+uorigin_p); |
642 687 | Int vcell=Int(vscale_p*v+vorigin_p); |
643 688 | imWeight(chn,row)=weight(chn%nChanWt,row); |
644 - | if((ucell>0)&&(ucell<nx_p)&&(vcell>0)&&(vcell<ny_p) &&gwt_p[fid](ucell,vcell)>0.0) { |
645 - | imWeight(chn,row)/=gwt_p[fid](ucell,vcell)*f2_p[fid]+d2_p[fid]; |
689 + | if((ucell>0)&&(ucell<nx_p)&&(vcell>0)&&(vcell<ny_p) &&a_gwt_p(ucell,vcell)>0.0) { |
690 + | imWeight(chn,row)/=a_gwt_p(ucell,vcell)*f2_p[fid]+d2_p[fid]; |
646 691 | sumwt+=imWeight(chn,row); |
647 692 | |
648 693 | } |
649 694 | else { |
650 695 | imWeight(chn,row)=0.0; |
651 696 | //ndrop++; |
652 697 | } |
653 698 | } |
654 699 | else{ |
655 700 | imWeight(chn,row)=0.0; |
727 772 | |
728 773 | |
729 774 | } |
730 775 | |
731 776 | VisImagingWeight& VisImagingWeight::operator=(const VisImagingWeight& other){ |
732 777 | if(this != &other){ |
733 778 | // gwt_p=other.gwt_p; |
734 779 | |
735 780 | gwt_p.resize(other.gwt_p.nelements(), true, false); |
736 781 | for (uInt k=0; k < gwt_p.nelements(); ++k){ |
737 - | gwt_p[k].resize(); |
782 + | //gwt_p[k].resize(); |
738 783 | gwt_p[k]=other.gwt_p[k]; |
739 784 | } |
740 785 | |
741 786 | |
742 787 | wgtType_p=other.wgtType_p; |
743 788 | uscale_p=other.uscale_p; |
744 789 | vscale_p=other.vscale_p; |
745 790 | f2_p.resize(); |
746 791 | d2_p.resize(); |
747 792 | f2_p=other.f2_p; |
767 812 | return wgtType_p; |
768 813 | |
769 814 | } |
770 815 | Bool VisImagingWeight::getWeightDensity (Block<Matrix<Float> >& density){ |
771 816 | if(wgtType_p != "uniform"){ |
772 817 | density.resize(0, true, false); |
773 818 | return false; |
774 819 | } |
775 820 | density.resize(gwt_p.nelements(), true, false); |
776 821 | for (uInt k=0; k < gwt_p.nelements(); ++k){ |
777 - | density[k].resize(); |
778 - | density[k].assign(gwt_p[k]); |
822 + | density[k].resize(gwt_p[k]->shape()); |
823 + | density[k]=gwt_p[k]->get(); |
779 824 | } |
780 825 | return true; |
781 826 | } |
782 827 | void VisImagingWeight::setWeightDensity(const Block<Matrix<Float> >& density){ |
783 828 | if(wgtType_p=="uniform"){ |
784 829 | gwt_p.resize(density.nelements(), true, false); |
785 830 | for (uInt k=0; k < gwt_p.nelements(); ++k){ |
786 - | gwt_p[k].resize(); |
787 - | gwt_p[k]=density[k]; |
831 + | //gwt_p[k].resize(); |
832 + | gwt_p[k]=new TempLattice<Float>(IPosition(2, density[k].shape()[0], density[k].shape()[1]),0); |
833 + | gwt_p[k]->put(density[k]); |
788 834 | } |
789 835 | //break any old reference |
790 836 | f2_p.resize(); |
791 837 | d2_p.resize(); |
792 838 | f2_p.resize(gwt_p.nelements()); |
793 839 | d2_p.resize(gwt_p.nelements()); |
794 840 | f2_p.set(1.0); |
795 841 | d2_p.set(0.0); |
796 842 | //Float f2, d2; |
797 843 | for(uInt fid=0; fid < gwt_p.nelements(); ++fid){ |
844 + | activeFieldIndex_p=fid; |
845 + | a_gwt_p=gwt_p[fid]->get(); |
798 846 | if (rmode_p=="norm") { |
799 847 | Double sumlocwt = 0.; |
800 - | for(Int vgrid=0;vgrid<gwt_p[fid].shape()(1);vgrid++) { |
801 - | for(Int ugrid=0;ugrid<gwt_p[fid].shape()(0);ugrid++) { |
802 - | if(gwt_p[fid](ugrid, vgrid)>0.0) sumlocwt+=square(gwt_p[fid](ugrid,vgrid)); |
848 + | for(Int vgrid=0;vgrid<a_gwt_p.shape()(1);vgrid++) { |
849 + | for(Int ugrid=0;ugrid<a_gwt_p.shape()(0);ugrid++) { |
850 + | if(a_gwt_p(ugrid, vgrid)>0.0) sumlocwt+=square(a_gwt_p(ugrid,vgrid)); |
803 851 | } |
804 852 | } |
805 - | Double sumwt_fid=sum(gwt_p[fid]); |
853 + | Double sumwt_fid=sum(a_gwt_p); |
806 854 | f2_p[fid] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumwt_fid); |
807 855 | d2_p[fid] = 1.0; |
808 856 | } |
809 857 | else if (rmode_p=="abs") { |
810 858 | f2_p[fid] = square(robust_p); |
811 859 | d2_p[fid] = 2.0 * square(noise_p.get("Jy").getValue()); |
812 860 | |
813 861 | } |
814 862 | else { |
815 863 | f2_p[fid] = 1.0; |