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         888 : FrequencySelection::addCorrelationSlices (const Vector <Vector <Slice> > & slices)
      17             : {
      18         888 :     correlationSlices_p = slices;
      19         888 : }
      20             : 
      21             : void
      22     1474191 : FrequencySelection::filterByWindow (Int windowId) const
      23             : {
      24     1474191 :     filterWindowId_p = windowId;
      25     1474191 : }
      26             : 
      27             : Int
      28     1810466 : FrequencySelection::filterWindow() const
      29             : {
      30     1810466 :     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           0 : }
      51             : 
      52             : Vector <Slice>
      53      743352 : FrequencySelection::getCorrelationSlices (Int polarizationId) const
      54             : {
      55      743352 :     if (correlationSlices_p.nelements() == 0){
      56      741965 :         return Vector<Slice> (); // empty vector implies all elements
      57             :     }
      58             : 
      59        1387 :     Assert (polarizationId >= 0 && polarizationId < (int) correlationSlices_p.nelements());
      60             : 
      61        1387 :     return correlationSlices_p [polarizationId];
      62             : }
      63             : 
      64             : Int
      65   563989361 : FrequencySelection::getFrameOfReference () const
      66             : {
      67   563989361 :     return referenceFrame_p;
      68             : }
      69             : 
      70             : void
      71        4397 : FrequencySelectionUsingChannels::add (Int spectralWindow, Int firstChannel,
      72             :                                       Int nChannels, Int increment)
      73             : {
      74        4397 :     Assert (spectralWindow >= 0);
      75        4397 :     Assert (firstChannel >= 0);
      76        4397 :     Assert (nChannels > 0);
      77        4397 :     Assert (increment != 0 || nChannels == 1);
      78             : 
      79        4397 :     elements_p.push_back (Element (spectralWindow, firstChannel, nChannels, increment));
      80        4397 : }
      81             : 
      82        9396 : FrequencySelectionUsingChannels::FrequencySelectionUsingChannels (const FrequencySelectionUsingChannels & other)
      83             : : FrequencySelection (other),
      84        9396 :   elements_p (other.elements_p),
      85        9396 :   filtered_p (other.filtered_p),
      86       18792 :   refinements_p (nullptr)
      87             : {
      88        9396 :     if (other.refinements_p){
      89           0 :         refinements_p.reset (new FrequencySelectionUsingFrame (* other.refinements_p.get()));
      90             :     }
      91        9396 : }
      92             : 
      93             : void
      94         811 : 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         811 :     MSSelection & msSelection = const_cast <MSSelection &> (msSelectionConst);
     104             :         // Needed because many msSelection methods are not const for some reason
     105             : 
     106         811 :     Matrix<Int> channelList = msSelection.getChanList();
     107             : 
     108        1114 :     for (Int row = 0; row < (Int) channelList.nrow(); row ++){
     109             : 
     110         303 :         Int nChannels = channelList (row, StopChannel) - channelList (row, FirstChannel);
     111         303 :         nChannels = nChannels / channelList (row, Step) + 1;
     112             : 
     113         303 :         add (channelList (row, SpectralWindowId),
     114         303 :              channelList (row, FirstChannel),
     115             :              nChannels,
     116         303 :              channelList (row, Step));
     117             :     }
     118             : 
     119             :     // Extract and add the correlation selection.
     120             : 
     121         811 :     Vector <Vector<Slice> > correlationSlices;
     122         811 :     msSelection.getCorrSlices (correlationSlices, ms);
     123             : 
     124         811 :     addCorrelationSlices (correlationSlices);
     125             : 
     126         811 : }
     127             : 
     128             : FrequencySelectionUsingChannels::const_iterator
     129       12513 : FrequencySelectionUsingChannels::begin () const
     130             : {
     131       12513 :     filtered_p.clear();
     132             : 
     133       12513 :     for (Elements::const_iterator i = elements_p.begin();
     134       21180 :          i != elements_p.end();
     135        8667 :          i ++){
     136             : 
     137        8667 :         if (filterWindow () < 0 || i->spectralWindow_p == filterWindow()){
     138        2325 :             filtered_p.push_back (* i);
     139             :         }
     140             :     }
     141             : 
     142       12513 :     return filtered_p.begin();
     143             : }
     144             : 
     145             : 
     146             : FrequencySelection *
     147        9396 : FrequencySelectionUsingChannels::clone () const
     148             : {
     149        9396 :     return new FrequencySelectionUsingChannels (* this);
     150             : }
     151             : 
     152             : Bool
     153      245848 : FrequencySelectionUsingChannels::empty () const
     154             : {
     155      245848 :     return elements_p.empty();
     156             : }
     157             : 
     158             : 
     159             : FrequencySelectionUsingChannels::const_iterator
     160       14838 : FrequencySelectionUsingChannels::end () const
     161             : {
     162       14838 :     return filtered_p.end();
     163             : }
     164             : 
     165             : set<int>
     166        7591 : FrequencySelectionUsingChannels::getSelectedWindows () const
     167             : {
     168        7591 :     set <int> result;
     169        7591 :     for (Elements::const_iterator i = elements_p.begin();
     170       11910 :          i != elements_p.end();
     171        4319 :          i ++){
     172             : 
     173        4319 :         result.insert (i->spectralWindow_p);
     174             : 
     175             :     }
     176             : 
     177        7591 :     return result;
     178           0 : }
     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           0 :             Slice s = spwcFetcher (spectralWindow,
     230             :                                    refiningElement.getBeginFrequency(),
     231             :                                    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       12513 : FrequencySelectionUsingChannels::refinementNeeded () const
     251             : {
     252       12513 :     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           0 : }
     307             : 
     308        4183 : FrequencySelectionUsingFrame::FrequencySelectionUsingFrame (MFrequency::Types frameOfReference)
     309        4183 : : FrequencySelection (frameOfReference)
     310        4183 : {}
     311             : 
     312             : void
     313        3267 : FrequencySelectionUsingFrame::add (Int spectralWindow, Double bottomFrequency,
     314             :                                    Double topFrequency)
     315             : {
     316        3267 :     elements_p.push_back (Element (spectralWindow, bottomFrequency, topFrequency));
     317        3267 : }
     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      730839 : FrequencySelectionUsingFrame::begin () const
     329             : {
     330      730839 :     filtered_p.clear();
     331             : 
     332      730839 :     for (Elements::const_iterator i = elements_p.begin();
     333     1627405 :          i != elements_p.end();
     334      896566 :          i ++){
     335             : 
     336      896566 :         if (filterWindow () < 0 || i->spectralWindow_p == filterWindow()){
     337      741805 :             filtered_p.push_back (* i);
     338             :         }
     339             :     }
     340             : 
     341      730839 :     return filtered_p.begin();
     342             : }
     343             : 
     344             : FrequencySelection *
     345        5051 : FrequencySelectionUsingFrame::clone () const
     346             : {
     347        5051 :     return new FrequencySelectionUsingFrame (* this);
     348             : }
     349             : 
     350             : Bool
     351        7750 : FrequencySelectionUsingFrame::empty () const
     352             : {
     353        7750 :     return elements_p.empty();
     354             : }
     355             : 
     356             : 
     357             : FrequencySelectionUsingFrame::const_iterator
     358     1472644 : FrequencySelectionUsingFrame::end () const
     359             : {
     360     1472644 :     return filtered_p.end();
     361             : }
     362             : 
     363             : set<int>
     364        2495 : FrequencySelectionUsingFrame::getSelectedWindows () const
     365             : {
     366        2495 :     set<int> result;
     367        2495 :     for (Elements::const_iterator i = elements_p.begin();
     368        5863 :          i != elements_p.end();
     369        3368 :          i ++){
     370             : 
     371        3368 :         result.insert (i->spectralWindow_p);
     372             :     }
     373             : 
     374        2495 :     return result;
     375           0 : }
     376             : 
     377             : Double
     378      741805 : FrequencySelectionUsingFrame::Element::getBeginFrequency () const
     379             : {
     380      741805 :     return beginFrequency_p;
     381             : }
     382             : 
     383             : Double
     384      741805 : FrequencySelectionUsingFrame::Element::getEndFrequency () const
     385             : {
     386      741805 :     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           0 : }
     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          17 :          Vector<Int> outspw;
     420          17 :          Vector<Int> outstart;
     421          17 :          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          17 :  }
     460       13408 : FrequencySelections::FrequencySelections ()
     461       13408 : : filterWindow_p (-1)
     462       13408 : {}
     463             : 
     464        5999 : FrequencySelections::FrequencySelections (const FrequencySelections & other)
     465             : {
     466        5999 :     for (Selections::const_iterator s = other.selections_p.begin();
     467       10377 :          s != other.selections_p.end(); s++){
     468        4378 :         selections_p.push_back ((* s)->clone());
     469             :     }
     470             : 
     471        5999 :     filterWindow_p = other.filterWindow_p;
     472        5999 :     selectedWindows_p = other.selectedWindows_p;
     473        5999 : }
     474             : 
     475       38814 : FrequencySelections::~FrequencySelections ()
     476             : {
     477       19407 :     for (Selections::const_iterator s = selections_p.begin();
     478       33854 :          s != selections_p.end(); s++){
     479       14447 :         delete (* s);
     480             :     }
     481       19407 : }
     482             : 
     483             : void
     484       10069 : FrequencySelections::add (const FrequencySelection & selection)
     485             : {
     486       10069 :     if (! selections_p.empty()){
     487         168 :         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       10069 :     selections_p.push_back (selection.clone());
     498       10069 :     Int msIndex = selections_p.size() - 1;
     499       10069 :     set<int> windows = selection.getSelectedWindows();
     500             : 
     501       17622 :     for (set<int>::const_iterator w = windows.begin(); w != windows.end(); w++){
     502        7553 :         selectedWindows_p.insert (make_pair (msIndex, * w));
     503             :     }
     504             : 
     505       10069 : }
     506             : 
     507             : FrequencySelections *
     508        5999 : FrequencySelections::clone () const
     509             : {
     510        5999 :     return new FrequencySelections (* this);
     511             : }
     512             : 
     513             : const FrequencySelection &
     514   266011753 : FrequencySelections::get (Int msIndex) const
     515             : {
     516   266011753 :     if (selections_p.empty()){
     517           0 :         return defaultSelection_p;
     518             :     }
     519             : 
     520   266011753 :     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   266011753 :     return * selections_p [msIndex];
     524             : }
     525             : 
     526             : 
     527             : Int
     528   296503350 : FrequencySelections::getFrameOfReference () const
     529             : {
     530   296503350 :     if (selections_p.empty()){
     531           0 :         return FrequencySelection::ByChannel;
     532             :     }
     533             :     else {
     534   296503350 :         return selections_p.front()->getFrameOfReference();
     535             :     }
     536             : }
     537             : 
     538             : Bool
     539      253598 : FrequencySelections::isSpectralWindowSelected (Int msIndex, Int spectralWindowId) const
     540             : {
     541             :     // Empty selections means everything is selected
     542             : 
     543      253598 :     if (selections_p.empty()){
     544           0 :         return true;
     545             :     }
     546             : 
     547      253598 :     Assert (msIndex >= 0 && msIndex < (int) selections_p.size() && selections_p [msIndex] != 0);
     548             : 
     549      253598 :     if (selections_p [msIndex]->empty()){
     550      240734 :         return true;
     551             :     }
     552             : 
     553             :     SelectedWindows::const_iterator swi =
     554       12864 :         selectedWindows_p.find (make_pair (msIndex, spectralWindowId));
     555             : 
     556       12864 :     return swi != selectedWindows_p.end();
     557             : }
     558             : 
     559             : 
     560             : Int
     561       11981 : FrequencySelections::size () const
     562             : {
     563       11981 :     return (Int) selections_p.size();
     564             : }
     565             : 
     566             : } // end namespace vi
     567             : 
     568             : using namespace casacore;
     569             : } // end namespace casa