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