Line data Source code
1 : //# BriggsCubeWeight.cc: Implementation for Briggs weighting for cubes
2 : //# Copyright (C) 2018-2021
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 3 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 : //# https://www.gnu.org/licenses/
16 : //#
17 : //# You should have received a copy of the GNU General Public License
18 : //# along with this library; if not, write to the Free Software Foundation,
19 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20 : //#
21 : //# Queries concerning CASA should be submitted at
22 : //# https://help.nrao.edu
23 : //#
24 : //# Postal address: CASA Project Manager
25 : //# National Radio Astronomy Observatory
26 : //# 520 Edgemont Road
27 : //# Charlottesville, VA 22903-2475 USA
28 : //#
29 : //#
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Arrays/Matrix.h>
32 : #include <casacore/casa/Arrays/Slice.h>
33 : #include <casacore/casa/Arrays/Vector.h>
34 : #include <casacore/casa/Containers/Record.h>
35 : #include <casacore/tables/DataMan/IncrementalStMan.h>
36 : #include <casacore/tables/DataMan/StManAipsIO.h>
37 : #include <casacore/tables/DataMan/StandardStMan.h>
38 : #include <casacore/tables/DataMan/TiledShapeStMan.h>
39 : #include <casacore/tables/Tables/ArrColDesc.h>
40 : #include <casacore/tables/Tables/ArrayColumn.h>
41 : #include <casacore/tables/Tables/ScaColDesc.h>
42 : #include <casacore/tables/Tables/ScalarColumn.h>
43 : #include <casacore/tables/Tables/TableUtil.h>
44 : #include <msvis/MSVis/MSUtil.h>
45 : #include <msvis/MSVis/VisBuffer2.h>
46 : #include <msvis/MSVis/VisImagingWeight.h>
47 : #include <msvis/MSVis/VisibilityIterator2.h>
48 : #include <synthesis/TransformMachines2/BriggsCubeWeightor.h>
49 : #include <synthesis/TransformMachines2/FTMachine.h>
50 :
51 : namespace casa { //# CASA namespace
52 : namespace refim { //# namespace refactor imaging
53 :
54 : using namespace casacore;
55 : using namespace casa;
56 : using namespace casacore;
57 : using namespace casa::refim;
58 : using namespace casacore;
59 : using namespace casa::vi;
60 :
61 0 : BriggsCubeWeightor::BriggsCubeWeightor()
62 : : grids_p(0), ft_p(0), f2_p(0), d2_p(0), uscale_p(0), vscale_p(0),
63 : uorigin_p(0), vorigin_p(0), nx_p(0), ny_p(0), rmode_p(""), noise_p(0.0),
64 : robust_p(2), superUniformBox_p(0), multiField_p(False),
65 : initialized_p(False), refFreq_p(-1.0),
66 : freqInterpMethod_p(InterpolateArray1D<Double, Complex>::nearestNeighbour),
67 0 : fracBW_p(0.0), wgtTab_p(nullptr), imWgtColName_p("") {
68 :
69 0 : multiFieldMap_p.clear();
70 0 : }
71 :
72 38 : BriggsCubeWeightor::BriggsCubeWeightor(
73 : const String &rmode, const Quantity &noise, const Double robust,
74 38 : const Double &fracBW, const Int superUniformBox, const Bool multiField)
75 : : grids_p(0), ft_p(0), f2_p(0), d2_p(0), uscale_p(0), vscale_p(0),
76 : uorigin_p(0), vorigin_p(0), nx_p(0), ny_p(0), initialized_p(False),
77 : refFreq_p(-1.0),
78 : freqInterpMethod_p(InterpolateArray1D<Double, Complex>::nearestNeighbour),
79 38 : wgtTab_p(nullptr), imWgtColName_p("") {
80 :
81 38 : rmode_p = rmode;
82 38 : noise_p = noise;
83 38 : robust_p = robust;
84 38 : superUniformBox_p = superUniformBox;
85 38 : multiField_p = multiField;
86 38 : multiFieldMap_p.clear();
87 38 : fracBW_p = fracBW;
88 38 : }
89 :
90 0 : BriggsCubeWeightor::BriggsCubeWeightor(
91 : vi::VisibilityIterator2 &vi, const String &rmode, const Quantity &noise,
92 : const Double robust, const ImageInterface<Complex> &templateimage,
93 : const RecordInterface &inrec, const Double &fracBW,
94 0 : const Int superUniformBox, const Bool multiField)
95 0 : : wgtTab_p(nullptr) {
96 0 : rmode_p = rmode;
97 0 : noise_p = noise;
98 0 : robust_p = robust;
99 0 : superUniformBox_p = superUniformBox;
100 0 : multiField_p = multiField;
101 0 : initialized_p = False;
102 0 : refFreq_p = -1.0;
103 0 : fracBW_p = fracBW;
104 :
105 0 : init(vi, templateimage, inrec);
106 0 : }
107 :
108 27 : String BriggsCubeWeightor::initImgWeightCol(
109 : vi::VisibilityIterator2 &vi, const ImageInterface<Complex> &templateimage,
110 : const RecordInterface &inRec) {
111 :
112 : // cout << "BriggsCubeWeightor::initImgWeightCol " << endl;
113 :
114 54 : CoordinateSystem cs = templateimage.coordinates();
115 :
116 54 : if (initialized_p && nx_p == templateimage.shape()(0) &&
117 27 : ny_p == templateimage.shape()(1)) {
118 :
119 0 : Double freq = cs.toWorld(IPosition(4, 0, 0, 0, 0))[3];
120 0 : if (freq == refFreq_p)
121 0 : return imWgtColName_p;
122 : }
123 :
124 27 : refFreq_p = cs.toWorld(IPosition(4, 0, 0, 0, 0))[3];
125 27 : visWgt_p = vi.getImagingWeightGenerator();
126 54 : VisImagingWeight vWghtNat("natural");
127 27 : vi.useImagingWeight(vWghtNat);
128 27 : vi::VisBuffer2 *vb = vi.getVisBuffer();
129 54 : std::vector<pair<Int, Int>> fieldsToUse;
130 54 : std::set<Int> msInUse;
131 27 : rownr_t nrows = 0;
132 54 : String ephemtab("");
133 27 : if (inRec.isDefined("ephemeristable")) {
134 0 : inRec.get("ephemeristable", ephemtab);
135 : }
136 27 : Double freqbeg = cs.toWorld(IPosition(4, 0, 0, 0, 0))[3];
137 : Double freqend =
138 27 : cs.toWorld(IPosition(4, 0, 0, 0, templateimage.shape()[3] - 1))[3];
139 :
140 80 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
141 12515 : for (vi.origin(); vi.more(); vi.next()) {
142 12462 : nrows += vb->nRows();
143 12462 : msInUse.insert(vb->msId());
144 12462 : if (multiField_p) {
145 :
146 2592 : pair<Int, Int> ms_field = make_pair(vb->msId(), vb->fieldId()[0]);
147 2592 : if (std::find(fieldsToUse.begin(), fieldsToUse.end(), ms_field) ==
148 5184 : fieldsToUse.end()) {
149 8 : fieldsToUse.push_back(ms_field);
150 : }
151 : }
152 : }
153 : }
154 27 : wgtTab_p = nullptr;
155 : Int tmpInt;
156 27 : inRec.get("freqinterpmethod", tmpInt);
157 27 : freqInterpMethod_p =
158 27 : static_cast<InterpolateArray1D<Double, Complex>::InterpolationMethod>(
159 : tmpInt);
160 54 : ostringstream oss;
161 27 : oss << std::setprecision(12) << nrows << "_" << freqbeg << "_" << freqend
162 27 : << "_" << rmode_p << "_" << robust_p << "_interp_" << tmpInt;
163 :
164 : // cerr << "STRING " << oss.str() << endl;;
165 27 : imWgtColName_p = makeScratchImagingWeightTable(wgtTab_p, oss.str());
166 27 : if (wgtTab_p->nrow() == nrows) {
167 : // cerr << "REUSING "<< wgtTab_p << endl;
168 16 : return imWgtColName_p;
169 : } else {
170 11 : wgtTab_p->lock(True, 100);
171 11 : wgtTab_p->removeRow(wgtTab_p->rowNumbers());
172 : }
173 : // cerr << "CREATING "<< wgtTab_p << endl;
174 11 : vbrowms2wgtrow_p.clear();
175 11 : if (fieldsToUse.size() == 0)
176 8 : fieldsToUse.push_back(make_pair(Int(-1), Int(-1)));
177 : // cerr << "FIELDs to use " << Vector<pair<Int,Int> >(fieldsToUse) << endl;
178 :
179 11 : Bool inOneGo = True;
180 11 : uInt allSwingPad = 4;
181 : // if multifield each field weight density are independent so no other field
182 : // or ms fields would contribute
183 11 : if (!multiField_p) {
184 8 : allSwingPad =
185 8 : estimateSwingChanPad(vi, -1, cs, templateimage.shape()[3], ephemtab);
186 : }
187 11 : if (msInUse.size() > 1) {
188 :
189 0 : if ((allSwingPad > 4) &&
190 0 : (allSwingPad > uInt(templateimage.shape()[3] / 10)))
191 0 : inOneGo = False;
192 : // cerr << "allSwingPad " << allSwingPad << " inOneGo " << inOneGo << endl;
193 : }
194 : ///////////////
195 : // cerr << "###fieldsInUSE " << Vector<pair<Int, Int> >(fieldsToUse) << endl;;
196 : // inOneGo=True;
197 :
198 : /////////////////////////////
199 11 : if (inOneGo) {
200 33 : fillImgWeightCol(vi, inRec, -1, fieldsToUse, allSwingPad,
201 22 : templateimage.shape(), cs);
202 : } else {
203 : /// Lets process the ms independently as swingpad can become very large for
204 : /// MSs seperated by large epochs
205 0 : for (auto msiter = msInUse.begin(); msiter != msInUse.end(); ++msiter) {
206 0 : uInt swingpad = estimateSwingChanPad(vi, *msiter, cs,
207 0 : templateimage.shape()[3], ephemtab);
208 0 : fillImgWeightCol(vi, inRec, *msiter, fieldsToUse, swingpad,
209 0 : templateimage.shape(), cs);
210 : }
211 : }
212 11 : initialized_p = True;
213 11 : wgtTab_p->unlock();
214 11 : return imWgtColName_p;
215 : }
216 :
217 11 : void BriggsCubeWeightor::fillImgWeightCol(
218 : vi::VisibilityIterator2 &vi, const Record &inRec, const Int msid,
219 : std::vector<pair<Int, Int>> &fieldsToUse, const uInt swingpad,
220 : const IPosition &origShp, CoordinateSystem csOrig) {
221 11 : vi::VisBuffer2 *vb = vi.getVisBuffer();
222 25 : for (uInt k = 0; k < fieldsToUse.size(); ++k) {
223 14 : vi.originChunks();
224 14 : vi.origin();
225 : // cerr << "####nchannels " << vb->nChannels() << " swingpad " << swingpad
226 : // << endl;
227 28 : IPosition shp = origShp;
228 14 : nx_p = shp[0];
229 14 : ny_p = shp[1];
230 28 : CoordinateSystem cs = csOrig;
231 : // CoordinateSystem cs=templateimage.coordinates();
232 :
233 28 : Vector<String> units = cs.worldAxisUnits();
234 14 : units[0] = "rad";
235 14 : units[1] = "rad";
236 14 : cs.setWorldAxisUnits(units);
237 28 : Vector<Double> incr = cs.increment();
238 14 : uscale_p = (nx_p * incr[0]);
239 14 : vscale_p = (ny_p * incr[1]);
240 14 : uorigin_p = nx_p / 2;
241 14 : vorigin_p = ny_p / 2;
242 : // IPosition shp=templateimage.shape();
243 14 : shp[3] = shp[3] + swingpad; // add extra channels at begining and end;
244 28 : Vector<Double> refpix = cs.referencePixel();
245 14 : refpix[3] += swingpad / 2;
246 14 : cs.setReferencePixel(refpix);
247 : // cerr << "REFPIX " << refpix << " new SHAPE " << shp << " swigpad " <<
248 : // swingpad << " msid "<< msid << " fieldToUse.msid "
249 : // <<fieldsToUse[k].first << " fieldid " << fieldsToUse[k].second << endl;
250 28 : TempImage<Complex> newTemplate(shp, cs);
251 28 : Matrix<Double> sumWgts;
252 :
253 14 : initializeFTMachine(0, newTemplate, inRec);
254 28 : Matrix<Float> dummy;
255 : // cerr << "new template shape " << newTemplate.shape() << endl;
256 14 : ft_p[0]->initializeToSky(newTemplate, dummy, *vb);
257 28 : Vector<Double> convFunc(2 + superUniformBox_p, 1.0);
258 : // cerr << "superuniform box " << superUniformBox_p << endl;
259 14 : ft_p[0]->modifyConvFunc(convFunc, superUniformBox_p, 1);
260 60 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
261 7384 : for (vi.origin(); vi.more(); vi.next()) {
262 : // cerr << "key and index "<< key << " " << index << " " <<
263 : // multiFieldMap_p[key] << endl;
264 7338 : if ((vb->fieldId()[0] == fieldsToUse[k].second &&
265 12732 : vb->msId() == fieldsToUse[k].first) ||
266 5394 : fieldsToUse[k].first == -1) {
267 5394 : ft_p[0]->put(*vb, -1, true, FTMachine::PSF);
268 : }
269 : }
270 : }
271 28 : Array<Float> griddedWeight;
272 14 : ft_p[0]->getGrid(griddedWeight);
273 : // cerr << " griddedWeight Shape " << griddedWeight.shape() << endl;
274 : // grids_p[index]->put(griddedWeight.reform(newTemplate.shape()));
275 14 : sumWgts = ft_p[0]->getSumWeights();
276 : // cerr << "sumweight " << sumWgts[index] << endl;
277 : // clear the ftmachine
278 14 : ft_p[0]->finalizeToSky();
279 :
280 14 : Int nchan = newTemplate.shape()(3);
281 : // cerr << "rmode " << rmode_p << endl;
282 :
283 233 : for (uInt chan = 0; chan < uInt(nchan); ++chan) {
284 438 : IPosition start(4, 0, 0, 0, chan);
285 438 : IPosition end(4, nx_p - 1, ny_p - 1, 0, chan);
286 : Matrix<Float> gwt(
287 657 : griddedWeight(start, end).reform(IPosition(2, nx_p, ny_p)));
288 390 : if ((rmode_p == "norm" || rmode_p == "bwtaper") &&
289 171 : (sumWgts(0, chan) >
290 : 0.0)) { // See CAS-13021 for bwtaper algorithm details
291 : // os << "Normal robustness, robust = " << robust << LogIO::POST;
292 129 : Double sumlocwt = 0.;
293 25045 : for (Int vgrid = 0; vgrid < ny_p; vgrid++) {
294 7509988 : for (Int ugrid = 0; ugrid < nx_p; ugrid++) {
295 7485072 : if (gwt(ugrid, vgrid) > 0.0)
296 130429 : sumlocwt += square(gwt(ugrid, vgrid));
297 : }
298 : }
299 258 : f2_p[0][chan] = square(5.0 * pow(10.0, Double(-robust_p))) /
300 129 : (sumlocwt / (2 * sumWgts(0, chan)));
301 129 : d2_p[0][chan] = 1.0;
302 :
303 90 : } else if (rmode_p == "abs") {
304 : // os << "Absolute robustness, robust = " << robust << ", noise = "
305 : // << noise.get("Jy").getValue() << "Jy" << LogIO::POST;
306 0 : f2_p[0][chan] = square(robust_p);
307 0 : d2_p[0][chan] = 2.0 * square(noise_p.get("Jy").getValue());
308 :
309 : } else {
310 90 : f2_p[0][chan] = 1.0;
311 90 : d2_p[0][chan] = 0.0;
312 : }
313 :
314 : } // chan
315 :
316 : // std::ofstream myfile;
317 : // myfile.open (wgtTab_p->tableName()+".txt");
318 14 : vi.originChunks();
319 14 : vi.origin();
320 60 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
321 7384 : for (vi.origin(); vi.more(); vi.next()) {
322 7338 : if ((msid < 0) || ((vb->msId()) == (msid))) {
323 7338 : if ((vb->fieldId()[0] == fieldsToUse[k].second &&
324 12732 : vb->msId() == fieldsToUse[k].first) ||
325 5394 : fieldsToUse[k].first == -1) {
326 10788 : Matrix<Float> imweight;
327 : /*///////////
328 : std::vector<Int> cmap=(ft_p[0]->channelMap(*vb)).tovector();
329 : int n=0;
330 : std::for_each(cmap.begin(), cmap.end(), [&n, &myfile](int
331 : &val){if(val > -1)myfile << " " << n; n++; }); myfile << "----msid
332 : " << vb->msId() << endl;
333 : /////////////////////////////*/
334 5394 : getWeightUniform(griddedWeight, imweight, *vb);
335 5394 : rownr_t nRows = vb->nRows();
336 : // Int nChans=vb->nChannels();
337 10788 : Vector<uInt> msId(nRows, uInt(vb->msId()));
338 10788 : RowNumbers rowidsnr = vb->rowIds();
339 10788 : Vector<uInt> rowids(rowidsnr.nelements());
340 5394 : convertArray(rowids, rowidsnr);
341 5394 : wgtTab_p->addRow(nRows, False);
342 5394 : rownr_t endrow = wgtTab_p->nrow() - 1;
343 5394 : rownr_t beginrow = endrow - nRows + 1;
344 : // Slicer sl(IPosition(2,beginrow,0), IPosition(2,endrow,nChans-1),
345 : // Slicer::endIsLast);
346 10788 : ArrayColumn<Float> col(*wgtTab_p, "IMAGING_WEIGHT");
347 : // cerr << "sl length " << sl.length() << " array shape " <<
348 : // fakeweight.shape() << " col length " << col.nrow() << endl;
349 : // col.putColumnRange(sl, fakeweight);
350 : // cerr << "nrows " << nRows << " imweight.shape " <<
351 : // imweight.shape() << endl;
352 1257168 : for (rownr_t row = 0; row < nRows; ++row)
353 1251774 : col.put(beginrow + row, imweight.column(row));
354 :
355 10788 : Slicer sl2(IPosition(1, endrow - nRows + 1), IPosition(1, endrow),
356 10788 : Slicer::endIsLast);
357 10788 : ScalarColumn<uInt> col2(*wgtTab_p, "MSID");
358 5394 : col2.putColumnRange(sl2, msId);
359 10788 : ScalarColumn<uInt> col3(*wgtTab_p, "ROWID");
360 :
361 5394 : col3.putColumnRange(sl2, rowids);
362 : }
363 : }
364 : }
365 : }
366 : // myfile.close();
367 : }
368 :
369 : ////Lets reset vi before returning
370 11 : vi.originChunks();
371 11 : vi.origin();
372 11 : }
373 :
374 8 : Int BriggsCubeWeightor::estimateSwingChanPad(vi::VisibilityIterator2 &vi,
375 : const Int msid,
376 : const CoordinateSystem &cs,
377 : const Int imNChan,
378 : const String &ephemtab) {
379 :
380 8 : vi.originChunks();
381 8 : vi.origin();
382 8 : vi::VisBuffer2 *vb = vi.getVisBuffer();
383 8 : Double freqbeg = cs.toWorld(IPosition(4, 0, 0, 0, 0))[3];
384 8 : Double freqend = cs.toWorld(IPosition(4, 0, 0, 0, imNChan - 1))[3];
385 8 : Double freqincr = fabs(cs.increment()[3]);
386 :
387 : SpectralCoordinate spCoord =
388 16 : cs.spectralCoordinate(cs.findCoordinate(Coordinate::SPECTRAL));
389 8 : MFrequency::Types freqframe = spCoord.frequencySystem(True);
390 : ////If using Undefined there is no Doppler correction to do
391 : // so no need of padding frequency
392 8 : if(freqframe == MFrequency::Undefined)
393 1 : return 0;
394 7 : uInt swingpad = 16;
395 7 : Double swingFreq = 0.0;
396 7 : Double minFreq = 1e99;
397 7 : Double maxFreq = 0.0;
398 14 : std::vector<Double> localminfreq, localmaxfreq, firstchanfreq;
399 7 : Int msID = -1;
400 7 : Int fieldID = -1;
401 7 : Int spwID = -1;
402 : ////TESTOO
403 : // int my_cpu_id;
404 : // MPI_Comm_rank(MPI_COMM_WORLD, &my_cpu_id);
405 : ///////////////////
406 : ///////
407 14 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
408 3367 : for (vi.origin(); vi.more(); vi.next()) {
409 : // process for required msid
410 3360 : if ((msid < 0) || (msid == vb->msId())) {
411 3360 : if ((msID != vb->msId()) || (fieldID != vb->fieldId()(0))) {
412 7 : msID = vb->msId();
413 7 : fieldID = vb->fieldId()(0);
414 7 : spwID = vb->spectralWindows()(0);
415 : Double localBeg, localEnd;
416 7 : Double localNchan = imNChan > 1 ? Double(imNChan - 1) : 1.0;
417 7 : Double localStep = abs(freqend - freqbeg) / localNchan;
418 7 : if (freqbeg < freqend) {
419 6 : localBeg = freqbeg;
420 6 : localEnd = freqend;
421 : } else {
422 1 : localBeg = freqend;
423 1 : localEnd = freqbeg;
424 : }
425 14 : Vector<Int> spw, start, nchan;
426 7 : if (ephemtab.size() != 0) {
427 0 : MSUtil::getSpwInSourceFreqRange(spw, start, nchan, vb->ms(),
428 : localBeg, localEnd, localStep,
429 : ephemtab, fieldID);
430 : } else {
431 7 : MSUtil::getSpwInFreqRange(spw, start, nchan, vb->ms(), localBeg,
432 : localEnd, localStep, freqframe, fieldID);
433 : }
434 14 : for (uInt spwk = 0; spwk < spw.nelements(); ++spwk) {
435 7 : if (spw[spwk] == spwID) {
436 7 : Vector<Double> mschanfreq = (vb->subtableColumns())
437 7 : .spectralWindow()
438 14 : .chanFreq()(spw[spwk]);
439 7 : if (mschanfreq[start[spwk] + nchan[spwk] - 1] >
440 7 : mschanfreq[start[spwk]]) {
441 6 : localminfreq.push_back(mschanfreq[start[spwk]]);
442 6 : localmaxfreq.push_back(
443 6 : mschanfreq[start[spwk] + nchan[spwk] - 1]);
444 : } else {
445 1 : localminfreq.push_back(
446 1 : mschanfreq[start[spwk] + nchan[spwk] - 1]);
447 1 : localmaxfreq.push_back(mschanfreq[start[spwk]]);
448 : }
449 :
450 7 : firstchanfreq.push_back(min(mschanfreq));
451 : // if(mschanfreq[start[spwk]+nchan[spwk]-1] <
452 : // localminfreq[localminfreq.size()-1])
453 : // localminfreq[localminfreq.size()-1]=mschanfreq[start[spwk]+nchan[spwk]-1];
454 7 : if (minFreq > localminfreq[localminfreq.size() - 1])
455 7 : minFreq = localminfreq[localminfreq.size() - 1];
456 7 : if (maxFreq < localmaxfreq[localmaxfreq.size() - 1])
457 7 : maxFreq = localmaxfreq[localmaxfreq.size() - 1];
458 : }
459 : }
460 : }
461 :
462 : } // input msid
463 : }
464 : }
465 :
466 7 : auto itf = firstchanfreq.begin();
467 7 : auto itmax = localmaxfreq.begin();
468 7 : Double firstchanshift = 0.0;
469 7 : Double minfirstchan = min(Vector<Double>(firstchanfreq));
470 14 : for (auto itmin = localminfreq.begin(); itmin != localminfreq.end();
471 7 : ++itmin) {
472 7 : if (swingFreq < abs(*itmin - minFreq))
473 0 : swingFreq = abs(*itmin - minFreq);
474 7 : if (swingFreq < abs(*itmax - maxFreq))
475 0 : swingFreq = abs(*itmax - maxFreq);
476 7 : if (firstchanshift < abs(*itf - minfirstchan))
477 0 : firstchanshift = abs(*itf - minfirstchan);
478 7 : itf++;
479 7 : itmax++;
480 : }
481 7 : Int extrapad = max(min(4, Int(imNChan / 10)), 1);
482 7 : swingpad =
483 7 : 2 * (Int(std::ceil((swingFreq + firstchanshift) / freqincr)) + extrapad);
484 : // cerr <<" swingfreq " << (swingFreq/freqincr) << " firstchanshift " <<
485 : // (firstchanshift/freqincr) << " SWINGPAD " << swingpad << endl;
486 : ////////////////
487 7 : return swingpad;
488 : }
489 :
490 27 : String BriggsCubeWeightor::makeScratchImagingWeightTable(
491 : CountedPtr<Table> &weightTable, const String &filetag) {
492 :
493 : // String wgtname=File::newUniqueName(".", "IMAGING_WEIGHT").absoluteName();
494 54 : String wgtname = Path("IMAGING_WEIGHT_" + filetag).absoluteName();
495 27 : if (Table::isReadable(wgtname)) {
496 16 : weightTable = new Table(wgtname, Table::Update);
497 16 : if (weightTable->nrow() > 0)
498 16 : return wgtname;
499 0 : weightTable = nullptr;
500 0 : TableUtil::deleteTable(wgtname, False);
501 : }
502 :
503 22 : TableDesc td;
504 11 : uInt cache_val = 32768;
505 11 : td.comment() = "Imaging_weight";
506 11 : td.addColumn(ScalarColumnDesc<uInt>("MSID"));
507 11 : td.addColumn(ScalarColumnDesc<uInt>("ROWID"));
508 11 : td.addColumn(ArrayColumnDesc<Float>("IMAGING_WEIGHT", 1));
509 :
510 11 : td.defineHypercolumn("TiledImagingWeight", 2,
511 22 : stringToVector("IMAGING_WEIGHT"));
512 22 : SetupNewTable newtab(wgtname, td, Table::New);
513 22 : IncrementalStMan incrStMan("MS_ID", cache_val);
514 11 : newtab.bindColumn("MSID", incrStMan);
515 22 : StandardStMan aipsStMan("ROWID", cache_val / 4);
516 11 : newtab.bindColumn("ROWID", aipsStMan);
517 33 : TiledShapeStMan tiledStMan("TiledImagingWeight", IPosition(2, 50, 500));
518 11 : newtab.bindColumn("IMAGING_WEIGHT", tiledStMan);
519 11 : weightTable = new Table(newtab);
520 : // weightTable->markForDelete();
521 11 : return wgtname;
522 : }
523 :
524 0 : void BriggsCubeWeightor::init(vi::VisibilityIterator2 &vi,
525 : const ImageInterface<Complex> &templateimage,
526 : const RecordInterface &inRec) {
527 : // cout << "BriggsCubeWeightor::init " << endl;
528 0 : LogIO os(LogOrigin("BriggsCubeWeightor", "constructor", WHERE));
529 :
530 : // freqInterpMethod_p=interpMethod;
531 : // freqFrameValid_p=freqFrameValid;
532 : // chanchunk may call the same object
533 0 : if (initialized_p && nx_p == templateimage.shape()(0) &&
534 0 : ny_p == templateimage.shape()(1)) {
535 0 : CoordinateSystem cs = templateimage.coordinates();
536 0 : Double freq = cs.toWorld(IPosition(4, 0, 0, 0, 0))[3];
537 0 : if (freq == refFreq_p)
538 0 : return;
539 : }
540 : // cerr << "in bgwt init " << endl;
541 : // Need to save previous wieght scheme of vi
542 0 : visWgt_p = vi.getImagingWeightGenerator();
543 0 : VisImagingWeight vWghtNat("natural");
544 0 : vi.useImagingWeight(vWghtNat);
545 0 : vi::VisBuffer2 *vb = vi.getVisBuffer();
546 0 : Int nIndices = 0;
547 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
548 0 : for (vi.origin(); vi.more(); vi.next()) {
549 0 : String key = String::toString(vb->msId()) + "_" +
550 0 : String::toString(vb->fieldId()(0));
551 0 : Int index = 0;
552 0 : if (multiField_p) {
553 : // find how many indices will be needed
554 0 : index = multiFieldMap_p.size();
555 0 : if (multiFieldMap_p.count(key) < 1)
556 0 : multiFieldMap_p[key] = index;
557 0 : nIndices = multiFieldMap_p.size();
558 : } else {
559 0 : multiFieldMap_p[key] = 0;
560 0 : nIndices = 1;
561 : }
562 : }
563 : }
564 : // cerr << "nindices " << nIndices << endl;
565 0 : vi.originChunks();
566 0 : vi.origin();
567 : String key =
568 0 : String::toString(vb->msId()) + "_" + String::toString(vb->fieldId()(0));
569 0 : IPosition shp = templateimage.shape();
570 0 : nx_p = shp[0];
571 0 : ny_p = shp[1];
572 0 : CoordinateSystem cs = templateimage.coordinates();
573 0 : refFreq_p = cs.toWorld(IPosition(4, 0, 0, 0, 0))[3];
574 0 : Vector<String> units = cs.worldAxisUnits();
575 0 : units[0] = "rad";
576 0 : units[1] = "rad";
577 0 : cs.setWorldAxisUnits(units);
578 0 : Vector<Double> incr = cs.increment();
579 0 : uscale_p = (nx_p * incr[0]);
580 0 : vscale_p = (ny_p * incr[1]);
581 0 : uorigin_p = nx_p / 2;
582 0 : vorigin_p = ny_p / 2;
583 : ////TESTOO
584 : // IPosition shp=templateimage.shape();
585 0 : shp[3] = shp[3] + 4; // add two channel at begining and end;
586 0 : Vector<Double> refpix = cs.referencePixel();
587 0 : refpix[3] += 2;
588 0 : cs.setReferencePixel(refpix);
589 0 : TempImage<Complex> newTemplate(shp, cs);
590 : ///////////////////////
591 : // ImageInterface<Complex>& newTemplate=const_cast<ImageInterface<Complex>&
592 : // >(templateimage);
593 0 : Vector<Matrix<Double>> sumWgts(nIndices);
594 :
595 0 : for (int index = 0; index < nIndices; ++index) {
596 0 : initializeFTMachine(index, newTemplate, inRec);
597 0 : Matrix<Float> dummy;
598 :
599 0 : ft_p[index]->initializeToSky(newTemplate, dummy, *vb);
600 0 : Vector<Double> convFunc(2 + superUniformBox_p, 1.0);
601 : // cerr << "superuniform box " << superUniformBox_p << endl;
602 0 : ft_p[index]->modifyConvFunc(convFunc, superUniformBox_p, 1);
603 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
604 0 : for (vi.origin(); vi.more(); vi.next()) {
605 :
606 0 : key = String::toString(vb->msId()) + "_" +
607 0 : String::toString(vb->fieldId()(0));
608 :
609 : // cerr << "key and index "<< key << " " << index << " " <<
610 : // multiFieldMap_p[key] << endl;
611 0 : if (multiFieldMap_p[key] == index) {
612 0 : ft_p[index]->put(*vb, -1, true, FTMachine::PSF);
613 : }
614 : }
615 : }
616 0 : Array<Float> griddedWeight;
617 0 : ft_p[index]->getGrid(griddedWeight);
618 : // cerr << index << " griddedWeight Shape " << griddedWeight.shape() <<
619 : // endl;
620 0 : grids_p[index]->put(griddedWeight.reform(newTemplate.shape()));
621 0 : sumWgts[index] = ft_p[index]->getSumWeights();
622 : // cerr << "sumweight " << sumWgts[index] << endl;
623 : // clear the ftmachine
624 0 : ft_p[index]->finalizeToSky();
625 : }
626 : ////Lets reset vi before returning
627 0 : vi.originChunks();
628 0 : vi.origin();
629 :
630 0 : Int nchan = newTemplate.shape()(3);
631 0 : for (uInt index = 0; index < f2_p.nelements(); ++index) {
632 : // cerr << "rmode " << rmode_p << endl;
633 :
634 0 : for (uInt chan = 0; chan < uInt(nchan); ++chan) {
635 0 : IPosition start(4, 0, 0, 0, chan);
636 0 : IPosition shape(4, nx_p, ny_p, 1, 1);
637 0 : Array<Float> arr;
638 0 : grids_p[index]->getSlice(arr, start, shape, True);
639 0 : Matrix<Float> gwt(arr);
640 0 : if ((rmode_p == "norm" || rmode_p == "bwtaper") &&
641 0 : (sumWgts[index](0, chan) >
642 : 0.0)) { // See CAS-13021 for bwtaper algorithm details
643 : // os << "Normal robustness, robust = " << robust << LogIO::POST;
644 0 : Double sumlocwt = 0.;
645 0 : for (Int vgrid = 0; vgrid < ny_p; vgrid++) {
646 0 : for (Int ugrid = 0; ugrid < nx_p; ugrid++) {
647 0 : if (gwt(ugrid, vgrid) > 0.0)
648 0 : sumlocwt += square(gwt(ugrid, vgrid));
649 : }
650 : }
651 0 : f2_p[index][chan] = square(5.0 * pow(10.0, Double(-robust_p))) /
652 0 : (sumlocwt / (2 * sumWgts[index](0, chan)));
653 0 : d2_p[index][chan] = 1.0;
654 :
655 0 : } else if (rmode_p == "abs") {
656 : // os << "Absolute robustness, robust = " << robust << ", noise = "
657 : // << noise.get("Jy").getValue() << "Jy" << LogIO::POST;
658 0 : f2_p[index][chan] = square(robust_p);
659 0 : d2_p[index][chan] = 2.0 * square(noise_p.get("Jy").getValue());
660 :
661 : } else {
662 0 : f2_p[index][chan] = 1.0;
663 0 : d2_p[index][chan] = 0.0;
664 : }
665 :
666 : } // chan
667 : }
668 0 : initialized_p = True;
669 : }
670 :
671 12462 : void BriggsCubeWeightor::weightUniform(Matrix<Float> &imweight,
672 : const vi::VisBuffer2 &vb) {
673 :
674 : // cout << "BriggsCubeWeightor::weightUniform" << endl;
675 :
676 12462 : if (!wgtTab_p.null())
677 12462 : return readWeightColumn(imweight, vb);
678 :
679 0 : if (multiFieldMap_p.size() == 0)
680 0 : throw(AipsError("BriggsCubeWeightor has not been initialized"));
681 : String key =
682 0 : String::toString(vb.msId()) + "_" + String::toString(vb.fieldId()(0));
683 0 : Int index = multiFieldMap_p[key];
684 0 : Vector<Int> chanMap = ft_p[0]->channelMap(vb);
685 : // cerr << "weightuniform chanmap " << chanMap << endl;
686 : /// No matching channels
687 0 : if (max(chanMap) == -1)
688 0 : return;
689 0 : Int nvischan = vb.nChannels();
690 0 : rownr_t nRow = vb.nRows();
691 0 : Matrix<Double> uvw = vb.uvw();
692 0 : imweight.resize(nvischan, nRow);
693 0 : imweight.set(0.0);
694 :
695 0 : Matrix<Float> weight;
696 0 : VisImagingWeight::unPolChanWeight(weight, vb.weightSpectrum());
697 0 : Matrix<Bool> flag;
698 0 : cube2Matrix(vb.flagCube(), flag);
699 :
700 0 : Int nChanWt = weight.shape()(0);
701 0 : Double sumwt = 0.0;
702 : Float u, v;
703 : Double fracBW, nCellsBW, uvDistanceFactor;
704 0 : IPosition pos(4, 0);
705 :
706 0 : fracBW = fracBW_p;
707 0 : if (rmode_p == "bwtaper") // See CAS-13021 for bwtaper algorithm details
708 : {
709 0 : if (fracBW == 0.0) {
710 : throw(AipsError(
711 0 : "BriggsCubeWeightor fractional bandwith is not a valid value, 0.0."));
712 : }
713 : // cout << "BriggsCubeWeightor::weightUniform fracBW " << fracBW << endl;
714 : }
715 :
716 0 : for (rownr_t row = 0; row < nRow; row++) {
717 0 : for (Int chn = 0; chn < nvischan; chn++) {
718 0 : if ((!flag(chn, row)) && (chanMap(chn) > -1)) {
719 0 : pos(3) = chanMap(chn);
720 0 : Float f = vb.getFrequency(0, chn) / C::c;
721 0 : u = -uvw(0, row) * f;
722 0 : v = -uvw(1, row) * f;
723 0 : Int ucell = Int(std::round(uscale_p * u + uorigin_p));
724 0 : Int vcell = Int(std::round(vscale_p * v + vorigin_p));
725 0 : pos(0) = ucell;
726 0 : pos(1) = vcell;
727 :
728 0 : imweight(chn, row) = 0.0;
729 0 : if ((ucell > 0) && (ucell < nx_p) && (vcell > 0) && (vcell < ny_p)) {
730 0 : Float gwt = grids_p[index]->getAt(pos);
731 0 : if (gwt > 0) {
732 0 : imweight(chn, row) = weight(chn % nChanWt, row);
733 :
734 0 : if (rmode_p ==
735 : "bwtaper") { // See CAS-13021 for bwtaper algorithm details
736 0 : nCellsBW = fracBW *
737 0 : sqrt(pow(uscale_p * u, 2.0) + pow(vscale_p * v, 2.0));
738 0 : uvDistanceFactor = nCellsBW + 0.5;
739 0 : if (uvDistanceFactor < 1.5)
740 0 : uvDistanceFactor = (4.0 - nCellsBW) / (4.0 - 2.0 * nCellsBW);
741 0 : imweight(chn, row) /=
742 0 : gwt * f2_p[index][pos[3]] / uvDistanceFactor +
743 0 : d2_p[index][pos[3]];
744 : } else {
745 0 : imweight(chn, row) /=
746 0 : gwt * f2_p[index][pos[3]] + d2_p[index][pos[3]];
747 : }
748 :
749 0 : sumwt += imweight(chn, row);
750 : }
751 : }
752 : } else {
753 0 : imweight(chn, row) = 0.0;
754 : }
755 : }
756 : }
757 :
758 0 : if (visWgt_p.doFilter()) {
759 0 : visWgt_p.filter(imweight, flag, uvw, vb.getFrequencies(0), imweight);
760 : }
761 : }
762 12462 : void BriggsCubeWeightor::readWeightColumn(
763 : casacore::Matrix<casacore::Float> &imweight, const vi::VisBuffer2 &vb) {
764 :
765 12462 : if (vbrowms2wgtrow_p.size() == 0) {
766 81 : Vector<uInt> msids = ScalarColumn<uInt>(*wgtTab_p, "MSID").getColumn();
767 81 : Vector<rownr_t> msrowid(ScalarColumn<uInt>(*wgtTab_p, "ROWID").nrow());
768 27 : convertArray(msrowid, ScalarColumn<uInt>(*wgtTab_p, "ROWID").getColumn());
769 3518829 : for (uInt k = 0; k < msids.nelements(); ++k) {
770 3518802 : vbrowms2wgtrow_p[make_pair(msids[k], msrowid[k])] = k;
771 : }
772 : // cerr << "Map size " << vbrowms2wgtrow_p.size() << " max size " <<
773 : // vbrowms2wgtrow_p.max_size() << endl;
774 : }
775 12462 : imweight.resize(vb.nChannels(), vb.nRows());
776 12462 : uInt msidnow = vb.msId();
777 24924 : RowNumbers rowIds = vb.rowIds();
778 24924 : ArrayColumn<Float> imwgtcol(*wgtTab_p, "IMAGING_WEIGHT");
779 3531264 : for (rownr_t k = 0; k < (vb.nRows()); ++k) {
780 3518802 : rownr_t tabrow = vbrowms2wgtrow_p[make_pair(msidnow, rowIds[k])];
781 : // cerr << imwgtcol.get(tabrow).shape() << " " <<
782 : // imweight.column(k).shape() << endl;
783 3518802 : imweight.column(k) = imwgtcol.get(tabrow);
784 : }
785 12462 : }
786 5394 : void BriggsCubeWeightor::getWeightUniform(const Array<Float> &wgtDensity,
787 : Matrix<Float> &imweight,
788 : const vi::VisBuffer2 &vb) {
789 : // cout << "BriggsCubeWeightor::getWeightUniform" << endl;
790 5394 : Vector<Int> chanMap = ft_p[0]->channelMap(vb);
791 : // cerr << "weightuniform chanmap " << chanMap << endl;
792 : ////
793 : // std::vector<Int> cmap=chanMap.tovector();
794 : // int n=0;
795 : // std::for_each(cmap.begin(), cmap.end(), [&n](int &val){if(val > -1)cerr <<
796 : // " " << n; n++; }); cerr << "----msid " << vb.msId() << endl;
797 5394 : Int nvischan = vb.nChannels();
798 5394 : rownr_t nRow = vb.nRows();
799 5394 : Matrix<Double> uvw = vb.uvw();
800 5394 : imweight.resize(nvischan, nRow);
801 5394 : imweight.set(0.0);
802 : /// No matching channels
803 5394 : if (max(chanMap) == -1)
804 0 : return;
805 10788 : Matrix<Float> weight;
806 5394 : VisImagingWeight::unPolChanWeight(weight, vb.weightSpectrum());
807 10788 : Matrix<Bool> flag;
808 5394 : cube2Matrix(vb.flagCube(), flag);
809 :
810 5394 : Int nChanWt = weight.shape()(0);
811 5394 : Double sumwt = 0.0;
812 : Float u, v;
813 : Double fracBW, nCellsBW, uvDistanceFactor;
814 10788 : IPosition pos(4, 0);
815 :
816 5394 : fracBW = fracBW_p;
817 5394 : if (rmode_p == "bwtaper") // See CAS-13021 for bwtaper algorithm details
818 : {
819 : // cout << "BriggsCubeWeightor::getWeightUniform bwtaper" << f2_p[0] <<
820 : // endl;
821 0 : if (fracBW == 0.0) {
822 : throw(AipsError(
823 0 : "BriggsCubeWeightor fractional bandwith is not a valid value, 0.0."));
824 : }
825 : // cout << "BriggsCubeWeightor::weightUniform fracBW " << fracBW << endl;
826 : }
827 :
828 1257168 : for (rownr_t row = 0; row < nRow; row++) {
829 22657968 : for (Int chn = 0; chn < nvischan; chn++) {
830 21406194 : if ((!flag(chn, row)) && (chanMap(chn) > -1)) {
831 20363409 : pos(3) = chanMap(chn);
832 20363409 : Float f = vb.getFrequency(0, chn) / C::c;
833 20363409 : u = -uvw(0, row) * f;
834 20363409 : v = -uvw(1, row) * f;
835 20363409 : Int ucell = Int(std::round(uscale_p * u + uorigin_p));
836 20363409 : Int vcell = Int(std::round(vscale_p * v + vorigin_p));
837 20363409 : pos(0) = ucell;
838 20363409 : pos(1) = vcell;
839 : ////TESTOO
840 : // if(row==0){
841 :
842 : // ofstream myfile;
843 : // myfile.open ("briggsLoc.txt", ios::out | ios::app | ios::ate );
844 : // myfile << vb.rowIds()(0) << " uv " << uvw.column(0) << " loc " <<
845 : // pos[0] << ", " << pos[1] << "\n"<< endl; myfile.close();
846 :
847 : //}
848 : //////
849 20363409 : imweight(chn, row) = 0.0;
850 20363409 : if ((ucell > 0) && (ucell < nx_p) && (vcell > 0) && (vcell < ny_p)) {
851 20363409 : Float gwt = wgtDensity(pos);
852 20363409 : if (gwt > 0) {
853 20227359 : imweight(chn, row) = weight(chn % nChanWt, row);
854 :
855 20227359 : if (rmode_p ==
856 : "bwtaper") { // See CAS-13021 for bwtaper algorithm details
857 0 : nCellsBW = fracBW *
858 0 : sqrt(pow(uscale_p * u, 2.0) + pow(vscale_p * v, 2.0));
859 0 : uvDistanceFactor = nCellsBW + 0.5;
860 0 : if (uvDistanceFactor < 1.5)
861 0 : uvDistanceFactor = (4.0 - nCellsBW) / (4.0 - 2.0 * nCellsBW);
862 0 : imweight(chn, row) /=
863 0 : gwt * f2_p[0][pos[3]] / uvDistanceFactor + d2_p[0][pos[3]];
864 : } else {
865 20227359 : imweight(chn, row) /= gwt * f2_p[0][pos[3]] + d2_p[0][pos[3]];
866 : }
867 :
868 20227359 : sumwt += imweight(chn, row);
869 : }
870 : }
871 : // else {
872 : // imweight(chn,row)=0.0;
873 : // ndrop++;
874 : //}
875 : } else {
876 1042785 : imweight(chn, row) = 0.0;
877 : }
878 : }
879 : }
880 :
881 5394 : if (visWgt_p.doFilter()) {
882 1440 : visWgt_p.filter(imweight, flag, uvw, vb.getFrequencies(0), imweight);
883 : }
884 : }
885 :
886 14 : void BriggsCubeWeightor::initializeFTMachine(
887 : const uInt index, const ImageInterface<Complex> &templateimage,
888 : const RecordInterface &inRec) {
889 14 : Int nchan = templateimage.shape()(3);
890 14 : if (ft_p.nelements() <= index) {
891 11 : ft_p.resize(index + 1);
892 11 : grids_p.resize(index + 1);
893 11 : f2_p.resize(index + 1);
894 11 : d2_p.resize(index + 1);
895 11 : f2_p[index] = Vector<Float>(nchan, 0.0);
896 11 : d2_p[index] = Vector<Float>(nchan, 0.0);
897 : }
898 14 : ft_p[index] =
899 28 : new refim::GridFT(Long(1000000), Int(200), "BOX", 1.0, true, false);
900 : Int tmpInt;
901 14 : inRec.get("freqinterpmethod", tmpInt);
902 14 : freqInterpMethod_p =
903 14 : static_cast<InterpolateArray1D<Double, Complex>::InterpolationMethod>(
904 : tmpInt);
905 14 : ft_p[index]->setFreqInterpolation(freqInterpMethod_p);
906 14 : inRec.get("freqframevalid", freqFrameValid_p);
907 14 : ft_p[index]->setFrameValidity(freqFrameValid_p);
908 14 : String error;
909 14 : if (!(ft_p[index]->recoverMovingSourceState(error, inRec)))
910 : throw(AipsError(
911 0 : "BriggsCubeWeightor could not get the state of the ftmachine:" +
912 0 : error));
913 : // remember to make the stokes I
914 42 : grids_p[index] = new TempImage<Float>(templateimage.shape(),
915 28 : templateimage.coordinates(), 0.0);
916 14 : }
917 5394 : void BriggsCubeWeightor::cube2Matrix(const Cube<Bool> &fcube,
918 : Matrix<Bool> &fMat) {
919 5394 : fMat.resize(fcube.shape()[1], fcube.shape()[2]);
920 : Bool deleteIt1;
921 : Bool deleteIt2;
922 5394 : const Bool *pcube = fcube.getStorage(deleteIt1);
923 5394 : Bool *pflags = fMat.getStorage(deleteIt2);
924 1257168 : for (uInt row = 0; row < fcube.shape()[2]; row++) {
925 22657968 : for (Int chn = 0; chn < fcube.shape()[1]; chn++) {
926 21406194 : *pflags = *pcube++;
927 42812388 : for (Int pol = 1; pol < fcube.shape()[0]; pol++, pcube++) {
928 21406194 : *pflags = *pcube ? *pcube : *pflags;
929 : }
930 21406194 : pflags++;
931 : }
932 : }
933 5394 : fcube.freeStorage(pcube, deleteIt1);
934 5394 : fMat.putStorage(pflags, deleteIt2);
935 5394 : }
936 :
937 : } // namespace refim
938 : } // namespace casa
|
