Line data Source code
1 : #include <stdcasa/UtilJ.h>
2 : #include <msvis/MSVis/ViFrequencySelection.h>
3 : #include <msvis/MSVis/MSUtil.h>
4 : #include <casacore/ms/MSSel/MSSelection.h>
5 :
6 : #include <utility>
7 :
8 : using namespace std;
9 :
10 : using namespace casacore;
11 : namespace casa { //# NAMESPACE CASA - BEGIN
12 :
13 : namespace vi {
14 :
15 : void
16 886 : FrequencySelection::addCorrelationSlices (const Vector <Vector <Slice> > & slices)
17 : {
18 886 : correlationSlices_p = slices;
19 886 : }
20 :
21 : void
22 1416679 : FrequencySelection::filterByWindow (Int windowId) const
23 : {
24 1416679 : filterWindowId_p = windowId;
25 1416679 : }
26 :
27 : Int
28 1720714 : FrequencySelection::filterWindow() const
29 : {
30 1720714 : return filterWindowId_p;
31 : }
32 :
33 : String
34 0 : FrequencySelection::frameName (Int referenceFrame)
35 : {
36 0 : String result;
37 :
38 0 : if (referenceFrame >= 0 && referenceFrame < MFrequency::N_Types){
39 :
40 0 : result = MFrequency::showType (referenceFrame);
41 : }
42 0 : else if (referenceFrame == ByChannel){
43 : }
44 : else{
45 :
46 0 : ThrowIf (true, String::format ("Unknown frame of reference: id=%d", referenceFrame));
47 : }
48 :
49 0 : return result;
50 : }
51 :
52 : Vector <Slice>
53 714568 : FrequencySelection::getCorrelationSlices (Int polarizationId) const
54 : {
55 714568 : if (correlationSlices_p.nelements() == 0){
56 713187 : return Vector<Slice> (); // empty vector implies all elements
57 : }
58 :
59 1381 : Assert (polarizationId >= 0 && polarizationId < (int) correlationSlices_p.nelements());
60 :
61 1381 : return correlationSlices_p [polarizationId];
62 : }
63 :
64 : Int
65 541716800 : FrequencySelection::getFrameOfReference () const
66 : {
67 541716800 : return referenceFrame_p;
68 : }
69 :
70 : void
71 4341 : FrequencySelectionUsingChannels::add (Int spectralWindow, Int firstChannel,
72 : Int nChannels, Int increment)
73 : {
74 4341 : Assert (spectralWindow >= 0);
75 4341 : Assert (firstChannel >= 0);
76 4341 : Assert (nChannels > 0);
77 4341 : Assert (increment != 0 || nChannels == 1);
78 :
79 4341 : elements_p.push_back (Element (spectralWindow, firstChannel, nChannels, increment));
80 4341 : }
81 :
82 9261 : FrequencySelectionUsingChannels::FrequencySelectionUsingChannels (const FrequencySelectionUsingChannels & other)
83 : : FrequencySelection (other),
84 9261 : elements_p (other.elements_p),
85 9261 : filtered_p (other.filtered_p),
86 9261 : refinements_p (nullptr)
87 : {
88 9261 : if (other.refinements_p){
89 0 : refinements_p.reset (new FrequencySelectionUsingFrame (* other.refinements_p.get()));
90 : }
91 9261 : }
92 :
93 : void
94 809 : FrequencySelectionUsingChannels::add (const MSSelection & msSelectionConst,
95 : const MeasurementSet * ms)
96 : {
97 : // Add in the frequency selection from the provided MSSelection object
98 : //
99 : // Meanings of columns in the "matrix" (actually used as a parallel array)
100 :
101 : enum {SpectralWindowId, FirstChannel, StopChannel, Step};
102 :
103 809 : MSSelection & msSelection = const_cast <MSSelection &> (msSelectionConst);
104 : // Needed because many msSelection methods are not const for some reason
105 :
106 1618 : Matrix<Int> channelList = msSelection.getChanList();
107 :
108 1109 : for (Int row = 0; row < (Int) channelList.nrow(); row ++){
109 :
110 300 : Int nChannels = channelList (row, StopChannel) - channelList (row, FirstChannel);
111 300 : nChannels = nChannels / channelList (row, Step) + 1;
112 :
113 300 : add (channelList (row, SpectralWindowId),
114 300 : channelList (row, FirstChannel),
115 : nChannels,
116 300 : channelList (row, Step));
117 : }
118 :
119 : // Extract and add the correlation selection.
120 :
121 1618 : Vector <Vector<Slice> > correlationSlices;
122 809 : msSelection.getCorrSlices (correlationSlices, ms);
123 :
124 809 : addCorrelationSlices (correlationSlices);
125 :
126 809 : }
127 :
128 : FrequencySelectionUsingChannels::const_iterator
129 12457 : FrequencySelectionUsingChannels::begin () const
130 : {
131 12457 : filtered_p.clear();
132 :
133 21014 : for (Elements::const_iterator i = elements_p.begin();
134 21014 : i != elements_p.end();
135 8557 : i ++){
136 :
137 8557 : if (filterWindow () < 0 || i->spectralWindow_p == filterWindow()){
138 2318 : filtered_p.push_back (* i);
139 : }
140 : }
141 :
142 12457 : return filtered_p.begin();
143 : }
144 :
145 :
146 : FrequencySelection *
147 9261 : FrequencySelectionUsingChannels::clone () const
148 : {
149 9261 : return new FrequencySelectionUsingChannels (* this);
150 : }
151 :
152 : Bool
153 245669 : FrequencySelectionUsingChannels::empty () const
154 : {
155 245669 : return elements_p.empty();
156 : }
157 :
158 :
159 : FrequencySelectionUsingChannels::const_iterator
160 14775 : FrequencySelectionUsingChannels::end () const
161 : {
162 14775 : return filtered_p.end();
163 : }
164 :
165 : set<int>
166 7461 : FrequencySelectionUsingChannels::getSelectedWindows () const
167 : {
168 7461 : set <int> result;
169 11724 : for (Elements::const_iterator i = elements_p.begin();
170 11724 : i != elements_p.end();
171 4263 : i ++){
172 :
173 4263 : result.insert (i->spectralWindow_p);
174 :
175 : }
176 :
177 7461 : return result;
178 : }
179 :
180 : Int
181 0 : FrequencySelectionUsingChannels::getNChannels (Int spectralWindowId) const
182 : {
183 :
184 0 : Int result = 0;
185 :
186 0 : if (elements_p.empty()){
187 :
188 : }
189 : else {
190 :
191 0 : for (Elements::const_iterator i = elements_p.begin();
192 0 : i != elements_p.end();
193 0 : i ++){
194 :
195 0 : if (i->spectralWindow_p == spectralWindowId){
196 0 : result += i->nChannels_p;
197 : }
198 : }
199 : }
200 0 : return result;
201 : }
202 :
203 : void
204 0 : FrequencySelectionUsingChannels::refine (const FrequencySelectionUsingFrame & refiningSelection)
205 : {
206 0 : refinements_p.reset (new FrequencySelectionUsingFrame (refiningSelection));
207 0 : }
208 :
209 : void
210 0 : FrequencySelectionUsingChannels::applyRefinement (std::function <casacore::Slice (int, double, double)> spwcFetcher) const
211 : {
212 : // Loop through each of the parts of the refining selection.
213 :
214 0 : for (auto refiningElement : * refinements_p.get()){
215 :
216 : // Get the spectral window and then look for parts of the original selection which
217 : // apply to the same spectral window.
218 :
219 0 : int spectralWindow = refiningElement.getSpectralWindow();
220 :
221 0 : for (auto & originalSelection : elements_p){
222 :
223 0 : if (originalSelection.spectralWindow_p != spectralWindow){
224 0 : continue; // wrong window, so skip it
225 : }
226 :
227 : // Convert the element into channels in the current spectral window.
228 :
229 : Slice s = spwcFetcher (spectralWindow,
230 : refiningElement.getBeginFrequency(),
231 0 : refiningElement.getEndFrequency());
232 :
233 0 : int firstRefiningChannel = s.start();
234 0 : int lastRefiningChannel = s.end();
235 :
236 : // Refine the original selection so that it only applies to the intersection between
237 : // the refining channel range and the original range.
238 : //
239 : // This only applies if the two intersect; otherwise the original selection element
240 : // is left unchanged.
241 :
242 0 : refineSelection (originalSelection, firstRefiningChannel, lastRefiningChannel);
243 : }
244 : }
245 :
246 0 : refinements_p.reset (nullptr); // All done so destroy them.
247 0 : }
248 :
249 : bool
250 12457 : FrequencySelectionUsingChannels::refinementNeeded () const
251 : {
252 12457 : return !! refinements_p;
253 : }
254 :
255 : void
256 0 : FrequencySelectionUsingChannels::refineSelection (FrequencySelectionUsingChannels::Element & originalElement,
257 : int firstRefiningChannel, int lastRefiningChannel) const
258 : {
259 0 : int originalLastChannel = originalElement.firstChannel_p +
260 0 : (originalElement.nChannels_p - 1) * originalElement.increment_p;
261 :
262 : // If the two ranges do not overlap, then skip this operation since there's not basis
263 : // for refinement. Presumably the refinement operation will apply to a different
264 : // selection range in this spectral window.
265 :
266 0 : if (firstRefiningChannel > originalLastChannel ||
267 0 : lastRefiningChannel < originalElement.firstChannel_p)
268 : {
269 0 : return; // No overlap, so refinement not possible
270 : }
271 :
272 : // Refine the channel interval of the existing selection.
273 :
274 0 : int newFirstChannel = max (originalElement.firstChannel_p, firstRefiningChannel);
275 0 : int newLastChannel = min (originalLastChannel, lastRefiningChannel);
276 :
277 0 : originalElement.firstChannel_p = newFirstChannel;
278 0 : originalElement.nChannels_p = (newLastChannel - newFirstChannel) / originalElement.increment_p + 1;
279 : }
280 :
281 : size_t
282 0 : FrequencySelectionUsingChannels::size () const
283 : {
284 0 : return elements_p.size();
285 : }
286 :
287 :
288 : String
289 0 : FrequencySelectionUsingChannels::toString () const
290 : {
291 0 : String s = String::format ("{frame='%s' {", frameName (getFrameOfReference()).c_str());
292 :
293 0 : for (Elements::const_iterator e = elements_p.begin();
294 0 : e != elements_p.end();
295 0 : e ++){
296 :
297 0 : s += String::format ("(spw=%d, 1st=%d, n=%d, inc=%d)",
298 0 : e->spectralWindow_p,
299 0 : e->firstChannel_p,
300 0 : e->nChannels_p,
301 0 : e->increment_p);
302 : }
303 0 : s += "}}";
304 :
305 0 : return s;
306 : }
307 :
308 3970 : FrequencySelectionUsingFrame::FrequencySelectionUsingFrame (MFrequency::Types frameOfReference)
309 3970 : : FrequencySelection (frameOfReference)
310 3970 : {}
311 :
312 : void
313 3041 : FrequencySelectionUsingFrame::add (Int spectralWindow, Double bottomFrequency,
314 : Double topFrequency)
315 : {
316 3041 : elements_p.push_back (Element (spectralWindow, bottomFrequency, topFrequency));
317 3041 : }
318 :
319 : //void
320 : //FrequencySelectionUsingFrame::add (Int spectralWindow, Double bottomFrequency,
321 : // Double topFrequency, Double increment)
322 : //{
323 : // elements_p.push_back (Elements (spectralWindow, bottomFrequency, topFrequency, increment));
324 : //}
325 :
326 :
327 : FrequencySelectionUsingFrame::const_iterator
328 702111 : FrequencySelectionUsingFrame::begin () const
329 : {
330 702111 : filtered_p.clear();
331 :
332 1553911 : for (Elements::const_iterator i = elements_p.begin();
333 1553911 : i != elements_p.end();
334 851800 : i ++){
335 :
336 851800 : if (filterWindow () < 0 || i->spectralWindow_p == filterWindow()){
337 712821 : filtered_p.push_back (* i);
338 : }
339 : }
340 :
341 702111 : return filtered_p.begin();
342 : }
343 :
344 : FrequencySelection *
345 4732 : FrequencySelectionUsingFrame::clone () const
346 : {
347 4732 : return new FrequencySelectionUsingFrame (* this);
348 : }
349 :
350 : Bool
351 7312 : FrequencySelectionUsingFrame::empty () const
352 : {
353 7312 : return elements_p.empty();
354 : }
355 :
356 :
357 : FrequencySelectionUsingFrame::const_iterator
358 1414932 : FrequencySelectionUsingFrame::end () const
359 : {
360 1414932 : return filtered_p.end();
361 : }
362 :
363 : set<int>
364 2343 : FrequencySelectionUsingFrame::getSelectedWindows () const
365 : {
366 2343 : set<int> result;
367 5470 : for (Elements::const_iterator i = elements_p.begin();
368 5470 : i != elements_p.end();
369 3127 : i ++){
370 :
371 3127 : result.insert (i->spectralWindow_p);
372 : }
373 :
374 2343 : return result;
375 : }
376 :
377 : Double
378 712821 : FrequencySelectionUsingFrame::Element::getBeginFrequency () const
379 : {
380 712821 : return beginFrequency_p;
381 : }
382 :
383 : Double
384 712821 : FrequencySelectionUsingFrame::Element::getEndFrequency () const
385 : {
386 712821 : return endFrequency_p;
387 : }
388 :
389 : int
390 0 : FrequencySelectionUsingFrame::Element::getSpectralWindow () const
391 : {
392 0 : return spectralWindow_p;
393 : }
394 :
395 : String
396 0 : FrequencySelectionUsingFrame::toString () const
397 : {
398 0 : String s = String::format ("{frame='%s' {", frameName (getFrameOfReference()).c_str());
399 :
400 0 : for (Elements::const_iterator e = elements_p.begin();
401 0 : e != elements_p.end();
402 0 : e ++){
403 :
404 0 : s += String::format ("(spw=%d, 1st=%g, n=%g, inc=%g)",
405 0 : e->spectralWindow_p,
406 0 : e->beginFrequency_p,
407 0 : e->endFrequency_p,
408 0 : e->increment_p);
409 : }
410 :
411 0 : s += "}}";
412 :
413 0 : return s;
414 : }
415 17 : std::map<int, std::pair<int, int> > FrequencySelectionUsingFrame::getChannelRange (const casacore::MeasurementSet& ms) const {
416 :
417 17 : map<int, pair<int, int> > retval;
418 17 : MFrequency::Types freqframe=MFrequency::castType(getFrameOfReference());
419 34 : Vector<Int> outspw;
420 34 : Vector<Int> outstart;
421 34 : Vector<Int> outnchan;
422 34 : for (Elements::const_iterator j = elements_p.begin(); j != elements_p.end(); j++){
423 : //cerr << "elements " << (j->spectralWindow_p) << " " << (j->beginFrequency_p)<< " " << (j->endFrequency_p)<< endl;
424 :
425 : Int outstart;
426 : Int outnchan;
427 17 : Int spw=j->spectralWindow_p;
428 17 : Double begFreq=j->beginFrequency_p;
429 17 : Double endFreq=j->endFrequency_p;
430 : //cerr << "spw " << spw << " BegFreq " << begFreq << " EndFreq " << endFreq << endl;
431 17 : MSUtil::getChannelRangeFromFreqRange(outstart,
432 : outnchan, ms, spw, begFreq, endFreq, 1.0e-6,freqframe);
433 17 : if(outstart > -1){
434 17 : auto old=retval.find(spw);
435 17 : if(old != retval.end()){
436 0 : pair<int, int> oldrange=old->second;
437 0 : if(oldrange.second >= outstart){
438 0 : if(outnchan > (oldrange.second-outstart+1+oldrange.first))
439 0 : oldrange.first=outnchan;
440 : else
441 0 : oldrange.first=oldrange.second-outstart+1+oldrange.first;
442 0 : oldrange.second=outstart;
443 : }
444 : else{
445 0 : if((oldrange.first -(outstart - oldrange.second)) <outnchan){
446 0 : oldrange.first=outnchan+outstart-oldrange.second+1;
447 : }
448 : }
449 0 : old->second=oldrange;
450 : }
451 : else{
452 :
453 17 : retval[int(spw)]=make_pair(int(outnchan), int(outstart));
454 : }
455 : }
456 : }
457 34 : return retval;
458 :
459 : }
460 13019 : FrequencySelections::FrequencySelections ()
461 13019 : : filterWindow_p (-1)
462 13019 : {}
463 :
464 5781 : FrequencySelections::FrequencySelections (const FrequencySelections & other)
465 : {
466 9987 : for (Selections::const_iterator s = other.selections_p.begin();
467 9987 : s != other.selections_p.end(); s++){
468 4206 : selections_p.push_back ((* s)->clone());
469 : }
470 :
471 5781 : filterWindow_p = other.filterWindow_p;
472 5781 : selectedWindows_p = other.selectedWindows_p;
473 5781 : }
474 :
475 18800 : FrequencySelections::~FrequencySelections ()
476 : {
477 32793 : for (Selections::const_iterator s = selections_p.begin();
478 46786 : s != selections_p.end(); s++){
479 13993 : delete (* s);
480 : }
481 18800 : }
482 :
483 : void
484 9787 : FrequencySelections::add (const FrequencySelection & selection)
485 : {
486 9787 : if (! selections_p.empty()){
487 138 : ThrowIf (getFrameOfReference() != selection.getFrameOfReference(),
488 : String::format ("Frequency selection #%d has incompatible frame of reference %d:%s "
489 : "(!= %d:%s)",
490 : selections_p.size() + 1,
491 : selection.getFrameOfReference(),
492 : FrequencySelection::frameName (selection.getFrameOfReference()).c_str(),
493 : getFrameOfReference(),
494 : FrequencySelection::frameName (getFrameOfReference()).c_str()));
495 : }
496 :
497 9787 : selections_p.push_back (selection.clone());
498 9787 : Int msIndex = selections_p.size() - 1;
499 19574 : set<int> windows = selection.getSelectedWindows();
500 :
501 17044 : for (set<int>::const_iterator w = windows.begin(); w != windows.end(); w++){
502 7257 : selectedWindows_p.insert (make_pair (msIndex, * w));
503 : }
504 :
505 9787 : }
506 :
507 : FrequencySelections *
508 5781 : FrequencySelections::clone () const
509 : {
510 5781 : return new FrequencySelections (* this);
511 : }
512 :
513 : const FrequencySelection &
514 256127129 : FrequencySelections::get (Int msIndex) const
515 : {
516 256127129 : if (selections_p.empty()){
517 0 : return defaultSelection_p;
518 : }
519 :
520 256127129 : ThrowIf (msIndex < 0 || msIndex >= (int) selections_p.size(),
521 : String::format ("MS index, %d, out of bounds [0,%d]", msIndex, selections_p.size() - 1));
522 :
523 256127129 : return * selections_p [msIndex];
524 : }
525 :
526 :
527 : Int
528 284172940 : FrequencySelections::getFrameOfReference () const
529 : {
530 284172940 : if (selections_p.empty()){
531 0 : return FrequencySelection::ByChannel;
532 : }
533 : else {
534 284172940 : return selections_p.front()->getFrameOfReference();
535 : }
536 : }
537 :
538 : Bool
539 252981 : FrequencySelections::isSpectralWindowSelected (Int msIndex, Int spectralWindowId) const
540 : {
541 : // Empty selections means everything is selected
542 :
543 252981 : if (selections_p.empty()){
544 0 : return true;
545 : }
546 :
547 252981 : Assert (msIndex >= 0 && msIndex < (int) selections_p.size() && selections_p [msIndex] != 0);
548 :
549 252981 : if (selections_p [msIndex]->empty()){
550 240607 : return true;
551 : }
552 :
553 : SelectedWindows::const_iterator swi =
554 12374 : selectedWindows_p.find (make_pair (msIndex, spectralWindowId));
555 :
556 12374 : return swi != selectedWindows_p.end();
557 : }
558 :
559 :
560 : Int
561 11474 : FrequencySelections::size () const
562 : {
563 11474 : return (Int) selections_p.size();
564 : }
565 :
566 : } // end namespace vi
567 :
568 : using namespace casacore;
569 : } // end namespace casa
|
