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           0 : FrequencySelection::addCorrelationSlices (const Vector <Vector <Slice> > & slices)
      17             : {
      18           0 :     correlationSlices_p = slices;
      19           0 : }
      20             : 
      21             : void
      22        2843 : FrequencySelection::filterByWindow (Int windowId) const
      23             : {
      24        2843 :     filterWindowId_p = windowId;
      25        2843 : }
      26             : 
      27             : Int
      28         174 : FrequencySelection::filterWindow() const
      29             : {
      30         174 :     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        2843 : FrequencySelection::getCorrelationSlices (Int polarizationId) const
      54             : {
      55        2843 :     if (correlationSlices_p.nelements() == 0){
      56        2843 :         return Vector<Slice> (); // empty vector implies all elements
      57             :     }
      58             : 
      59           0 :     Assert (polarizationId >= 0 && polarizationId < (int) correlationSlices_p.nelements());
      60             : 
      61           0 :     return correlationSlices_p [polarizationId];
      62             : }
      63             : 
      64             : Int
      65      595148 : FrequencySelection::getFrameOfReference () const
      66             : {
      67      595148 :     return referenceFrame_p;
      68             : }
      69             : 
      70             : void
      71          14 : FrequencySelectionUsingChannels::add (Int spectralWindow, Int firstChannel,
      72             :                                       Int nChannels, Int increment)
      73             : {
      74          14 :     Assert (spectralWindow >= 0);
      75          14 :     Assert (firstChannel >= 0);
      76          14 :     Assert (nChannels > 0);
      77          14 :     Assert (increment != 0 || nChannels == 1);
      78             : 
      79          14 :     elements_p.push_back (Element (spectralWindow, firstChannel, nChannels, increment));
      80          14 : }
      81             : 
      82         849 : FrequencySelectionUsingChannels::FrequencySelectionUsingChannels (const FrequencySelectionUsingChannels & other)
      83             : : FrequencySelection (other),
      84         849 :   elements_p (other.elements_p),
      85         849 :   filtered_p (other.filtered_p),
      86        1698 :   refinements_p (nullptr)
      87             : {
      88         849 :     if (other.refinements_p){
      89           0 :         refinements_p.reset (new FrequencySelectionUsingFrame (* other.refinements_p.get()));
      90             :     }
      91         849 : }
      92             : 
      93             : void
      94           0 : 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           0 :     MSSelection & msSelection = const_cast <MSSelection &> (msSelectionConst);
     104             :         // Needed because many msSelection methods are not const for some reason
     105             : 
     106           0 :     Matrix<Int> channelList = msSelection.getChanList();
     107             : 
     108           0 :     for (Int row = 0; row < (Int) channelList.nrow(); row ++){
     109             : 
     110           0 :         Int nChannels = channelList (row, StopChannel) - channelList (row, FirstChannel);
     111           0 :         nChannels = nChannels / channelList (row, Step) + 1;
     112             : 
     113           0 :         add (channelList (row, SpectralWindowId),
     114           0 :              channelList (row, FirstChannel),
     115             :              nChannels,
     116           0 :              channelList (row, Step));
     117             :     }
     118             : 
     119             :     // Extract and add the correlation selection.
     120             : 
     121           0 :     Vector <Vector<Slice> > correlationSlices;
     122           0 :     msSelection.getCorrSlices (correlationSlices, ms);
     123             : 
     124           0 :     addCorrelationSlices (correlationSlices);
     125             : 
     126           0 : }
     127             : 
     128             : FrequencySelectionUsingChannels::const_iterator
     129        2843 : FrequencySelectionUsingChannels::begin () const
     130             : {
     131        2843 :     filtered_p.clear();
     132             : 
     133        2843 :     for (Elements::const_iterator i = elements_p.begin();
     134        2930 :          i != elements_p.end();
     135          87 :          i ++){
     136             : 
     137          87 :         if (filterWindow () < 0 || i->spectralWindow_p == filterWindow()){
     138          14 :             filtered_p.push_back (* i);
     139             :         }
     140             :     }
     141             : 
     142        2843 :     return filtered_p.begin();
     143             : }
     144             : 
     145             : 
     146             : FrequencySelection *
     147         849 : FrequencySelectionUsingChannels::clone () const
     148             : {
     149         849 :     return new FrequencySelectionUsingChannels (* this);
     150             : }
     151             : 
     152             : Bool
     153        6634 : FrequencySelectionUsingChannels::empty () const
     154             : {
     155        6634 :     return elements_p.empty();
     156             : }
     157             : 
     158             : 
     159             : FrequencySelectionUsingChannels::const_iterator
     160        2857 : FrequencySelectionUsingChannels::end () const
     161             : {
     162        2857 :     return filtered_p.end();
     163             : }
     164             : 
     165             : set<int>
     166         825 : FrequencySelectionUsingChannels::getSelectedWindows () const
     167             : {
     168         825 :     set <int> result;
     169         825 :     for (Elements::const_iterator i = elements_p.begin();
     170         839 :          i != elements_p.end();
     171          14 :          i ++){
     172             : 
     173          14 :         result.insert (i->spectralWindow_p);
     174             : 
     175             :     }
     176             : 
     177         825 :     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        2843 : FrequencySelectionUsingChannels::refinementNeeded () const
     251             : {
     252        2843 :     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           0 : FrequencySelectionUsingFrame::FrequencySelectionUsingFrame (MFrequency::Types frameOfReference)
     309           0 : : FrequencySelection (frameOfReference)
     310           0 : {}
     311             : 
     312             : void
     313           0 : FrequencySelectionUsingFrame::add (Int spectralWindow, Double bottomFrequency,
     314             :                                    Double topFrequency)
     315             : {
     316           0 :     elements_p.push_back (Element (spectralWindow, bottomFrequency, topFrequency));
     317           0 : }
     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           0 : FrequencySelectionUsingFrame::begin () const
     329             : {
     330           0 :     filtered_p.clear();
     331             : 
     332           0 :     for (Elements::const_iterator i = elements_p.begin();
     333           0 :          i != elements_p.end();
     334           0 :          i ++){
     335             : 
     336           0 :         if (filterWindow () < 0 || i->spectralWindow_p == filterWindow()){
     337           0 :             filtered_p.push_back (* i);
     338             :         }
     339             :     }
     340             : 
     341           0 :     return filtered_p.begin();
     342             : }
     343             : 
     344             : FrequencySelection *
     345           0 : FrequencySelectionUsingFrame::clone () const
     346             : {
     347           0 :     return new FrequencySelectionUsingFrame (* this);
     348             : }
     349             : 
     350             : Bool
     351           0 : FrequencySelectionUsingFrame::empty () const
     352             : {
     353           0 :     return elements_p.empty();
     354             : }
     355             : 
     356             : 
     357             : FrequencySelectionUsingFrame::const_iterator
     358           0 : FrequencySelectionUsingFrame::end () const
     359             : {
     360           0 :     return filtered_p.end();
     361             : }
     362             : 
     363             : set<int>
     364           0 : FrequencySelectionUsingFrame::getSelectedWindows () const
     365             : {
     366           0 :     set<int> result;
     367           0 :     for (Elements::const_iterator i = elements_p.begin();
     368           0 :          i != elements_p.end();
     369           0 :          i ++){
     370             : 
     371           0 :         result.insert (i->spectralWindow_p);
     372             :     }
     373             : 
     374           0 :     return result;
     375           0 : }
     376             : 
     377             : Double
     378           0 : FrequencySelectionUsingFrame::Element::getBeginFrequency () const
     379             : {
     380           0 :     return beginFrequency_p;
     381             : }
     382             : 
     383             : Double
     384           0 : FrequencySelectionUsingFrame::Element::getEndFrequency () const
     385             : {
     386           0 :     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           0 :  std::map<int, std::pair<int, int>  >  FrequencySelectionUsingFrame::getChannelRange (const casacore::MeasurementSet& ms) const {
     416             :          
     417           0 :          map<int, pair<int, int> > retval;
     418           0 :          MFrequency::Types freqframe=MFrequency::castType(getFrameOfReference());
     419           0 :          Vector<Int> outspw;
     420           0 :          Vector<Int> outstart;
     421           0 :          Vector<Int> outnchan;
     422           0 :          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           0 :                 Int spw=j->spectralWindow_p;
     428           0 :                 Double begFreq=j->beginFrequency_p;
     429           0 :                 Double endFreq=j->endFrequency_p;
     430             :                 //cerr << "spw " << spw << " BegFreq " << begFreq << "  EndFreq " << endFreq << endl;
     431           0 :                 MSUtil::getChannelRangeFromFreqRange(outstart,
     432             :                           outnchan, ms, spw, begFreq, endFreq, 1.0e-6,freqframe); 
     433           0 :                 if(outstart > -1){
     434           0 :                         auto  old=retval.find(spw);
     435           0 :                         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           0 :                                 retval[int(spw)]=make_pair(int(outnchan), int(outstart));
     454             :                         }
     455             :                 }
     456             :          }
     457           0 :          return retval;
     458             :          
     459           0 :  }
     460         825 : FrequencySelections::FrequencySelections ()
     461         825 : : filterWindow_p (-1)
     462         825 : {}
     463             : 
     464          24 : FrequencySelections::FrequencySelections (const FrequencySelections & other)
     465             : {
     466          24 :     for (Selections::const_iterator s = other.selections_p.begin();
     467          48 :          s != other.selections_p.end(); s++){
     468          24 :         selections_p.push_back ((* s)->clone());
     469             :     }
     470             : 
     471          24 :     filterWindow_p = other.filterWindow_p;
     472          24 :     selectedWindows_p = other.selectedWindows_p;
     473          24 : }
     474             : 
     475        1698 : FrequencySelections::~FrequencySelections ()
     476             : {
     477         849 :     for (Selections::const_iterator s = selections_p.begin();
     478        1698 :          s != selections_p.end(); s++){
     479         849 :         delete (* s);
     480             :     }
     481         849 : }
     482             : 
     483             : void
     484         825 : FrequencySelections::add (const FrequencySelection & selection)
     485             : {
     486         825 :     if (! selections_p.empty()){
     487           0 :         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         825 :     selections_p.push_back (selection.clone());
     498         825 :     Int msIndex = selections_p.size() - 1;
     499         825 :     set<int> windows = selection.getSelectedWindows();
     500             : 
     501         839 :     for (set<int>::const_iterator w = windows.begin(); w != windows.end(); w++){
     502          14 :         selectedWindows_p.insert (make_pair (msIndex, * w));
     503             :     }
     504             : 
     505         825 : }
     506             : 
     507             : FrequencySelections *
     508          24 : FrequencySelections::clone () const
     509             : {
     510          24 :     return new FrequencySelections (* this);
     511             : }
     512             : 
     513             : const FrequencySelection &
     514      527028 : FrequencySelections::get (Int msIndex) const
     515             : {
     516      527028 :     if (selections_p.empty()){
     517           0 :         return defaultSelection_p;
     518             :     }
     519             : 
     520      527028 :     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      527028 :     return * selections_p [msIndex];
     524             : }
     525             : 
     526             : 
     527             : Int
     528       65277 : FrequencySelections::getFrameOfReference () const
     529             : {
     530       65277 :     if (selections_p.empty()){
     531           0 :         return FrequencySelection::ByChannel;
     532             :     }
     533             :     else {
     534       65277 :         return selections_p.front()->getFrameOfReference();
     535             :     }
     536             : }
     537             : 
     538             : Bool
     539        6634 : FrequencySelections::isSpectralWindowSelected (Int msIndex, Int spectralWindowId) const
     540             : {
     541             :     // Empty selections means everything is selected
     542             : 
     543        6634 :     if (selections_p.empty()){
     544           0 :         return true;
     545             :     }
     546             : 
     547        6634 :     Assert (msIndex >= 0 && msIndex < (int) selections_p.size() && selections_p [msIndex] != 0);
     548             : 
     549        6634 :     if (selections_p [msIndex]->empty()){
     550        6604 :         return true;
     551             :     }
     552             : 
     553             :     SelectedWindows::const_iterator swi =
     554          30 :         selectedWindows_p.find (make_pair (msIndex, spectralWindowId));
     555             : 
     556          30 :     return swi != selectedWindows_p.end();
     557             : }
     558             : 
     559             : 
     560             : Int
     561           3 : FrequencySelections::size () const
     562             : {
     563           3 :     return (Int) selections_p.size();
     564             : }
     565             : 
     566             : } // end namespace vi
     567             : 
     568             : using namespace casacore;
     569             : } // end namespace casa