Changeset 344
- Timestamp:
- 02/01/05 16:17:25 (20 years ago)
- Location:
- trunk/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/SDLineFinder.cc
r343 r344 44 44 using namespace boost::python; 45 45 46 namespace asap { 47 48 /////////////////////////////////////////////////////////////////////////////// 49 // 50 // IStatHolder - an abstract class to collect statistics from the running 51 // mean calculator, if necessary. 52 // We define it here, because it is used in LFRunningMean and 53 // SDLineFinder only 54 // 55 56 struct IStatHolder { 57 // This function is called for each spectral channel processed by 58 // the running mean calculator. The order of channel numbers may be 59 // arbitrary 60 // ch - a number of channel, corresponding to (approximately) the centre 61 // of the running box 62 // box_nchan - number of channels in the box 63 // 64 virtual void accumulate(int ch, Float sum, Float sum2, int box_nchan) 65 throw(AipsError) = 0; 66 }; 67 68 // 69 /////////////////////////////////////////////////////////////////////////////// 70 71 /////////////////////////////////////////////////////////////////////////////// 72 // 73 // SignAccumulator - a simple class to deal with running mean statistics: 74 // it stores the sign of the value-mean only 75 // 76 class SignAccumulator : public IStatHolder { 77 Vector<Int> sign; // either +1, -1 or 0 78 const Vector<Float> &spectrum; // a reference to the spectrum 79 // to calculate the sign 80 public: 81 // all channels >=nchan are ignored 82 SignAccumulator(uInt nchan, const Vector<Float> &in_spectrum); 83 84 // This function is called for each spectral channel processed by 85 // the running mean calculator. The order of channel numbers may be 86 // arbitrary 87 // ch - a number of channel, corresponding to (approximately) the centre 88 // of the running box 89 // box_nchan - number of channels in the box 90 // 91 virtual void accumulate(int ch, Float sum, Float, int box_nchan) 92 throw(AipsError); 93 94 // access to the sign 95 const Vector<Int>& getSigns() const throw(); 96 }; 97 98 // 99 /////////////////////////////////////////////////////////////////////////////// 100 46 101 /////////////////////////////////////////////////////////////////////////////// 47 102 // … … 49 104 // 50 105 // 51 52 namespace asap {53 106 54 107 // An auxiliary class implementing one pass of the line search algorithm, … … 96 149 97 150 // find spectral lines and add them into list 98 void findLines(std::list<pair<int,int> > &lines) throw(casa::AipsError); 151 // if statholder is not NULL, the accumulate function of it will be 152 // called for each channel to save statistics 153 void findLines(std::list<pair<int,int> > &lines, 154 IStatHolder* statholder = NULL) throw(casa::AipsError); 99 155 100 156 protected: … … 123 179 124 180 }; 181 182 // 183 /////////////////////////////////////////////////////////////////////////////// 125 184 } // namespace asap 185 186 /////////////////////////////////////////////////////////////////////////////// 187 // 188 // SignAccumulator - a simple class to deal with running mean statistics: 189 // it stores the sign of the value-mean only 190 // 191 192 // all channels >=nchan are ignored 193 SignAccumulator::SignAccumulator(uInt nchan, 194 const Vector<Float> &in_spectrum) : sign(nchan,0), 195 spectrum(in_spectrum) {} 196 197 198 // This function is called for each spectral channel processed by 199 // the running mean calculator. The order of channel numbers may be 200 // arbitrary 201 // ch - a number of channel, corresponding to (approximately) the centre 202 // of the running box 203 // box_nchan - number of channels in the box 204 // 205 void SignAccumulator::accumulate(int ch, Float sum, Float, int box_nchan) 206 throw(AipsError) 207 { 208 if (ch>=sign.nelements()) return; 209 DebugAssert(ch>=0,AipsError); 210 DebugAssert(ch<=spectrum.nelements(), AipsError); 211 if (box_nchan) { 212 Float buf=spectrum[ch]-sum/Float(box_nchan); 213 if (buf>0) sign[ch]=1; 214 else if (buf<0) sign[ch]=-1; 215 else sign[ch]=0; 216 } else sign[ch]=0; 217 } 218 219 // access to the sign 220 const Vector<Int>& SignAccumulator::getSigns() const throw() 221 { 222 return sign; 223 } 224 126 225 127 226 // … … 253 352 254 353 // find spectral lines and add them into list 255 void LFRunningMean::findLines(std::list<pair<int,int> > &lines) 354 void LFRunningMean::findLines(std::list<pair<int,int> > &lines, 355 IStatHolder* statholder) 256 356 throw(casa::AipsError) 257 357 { … … 273 373 is_detected_before=False; 274 374 375 if (statholder!=NULL) 376 for (int n=0;n<initial_box_ch-max_box_nchan/2;++n) 377 statholder->accumulate(n,sum,sumsq,box_chan_cntr); 378 275 379 if (box_chan_cntr>=minboxnchan) 276 380 // there is a minimum amount of data. We can search in the … … 283 387 for (int n=initial_box_ch-max_box_nchan/2;n<edge.second;++n) { 284 388 advanceRunningBox(n+max_box_nchan/2); // update running mean & variance 389 if (statholder!=NULL) 390 statholder->accumulate(n,sum,sumsq,box_chan_cntr); 285 391 if (box_chan_cntr>=minboxnchan) // have enough data to process 286 392 processChannel(lines,n); 287 393 else processCurLine(lines); // just finish what was accumulated before 288 } 394 } 395 if (statholder!=NULL) 396 for (int n=edge.second;n<spectrum.nelements();++n) 397 statholder->accumulate(n,sum,sumsq,box_chan_cntr); 289 398 } 290 399 … … 468 577 lines.resize(0); // search from the scratch 469 578 Vector<Bool> temp_mask(mask); 470 Bool need2iterate=True;471 while ( need2iterate) {579 580 while (true) { 472 581 // line find algorithm 473 582 LFRunningMean lfalg(spectrum,temp_mask,edge,max_box_nchan,min_nchan, 474 583 threshold); 584 SignAccumulator sacc(spectrum.nelements(),spectrum); 585 586 // a buffer for new lines found at this iteration 475 587 std::list<pair<int,int> > new_lines; 476 lfalg.findLines(new_lines); 477 if (!new_lines.size()) need2iterate=False; 478 temp_mask=getMask(); 588 589 lfalg.findLines(new_lines,&sacc); 590 if (!new_lines.size()) break; // nothing new 591 592 searchForWings(new_lines, sacc.getSigns()); 593 479 594 // update the list (lines) merging intervals, if necessary 480 addNewSearchResult(new_lines); 595 addNewSearchResult(new_lines,lines); 596 // get a new mask 597 temp_mask=getMask(); 481 598 } 482 599 return int(lines.size()); … … 569 686 // concatenate two lists preserving the order. If two lines appear to 570 687 // be adjacent, they are joined into the new one 571 void SDLineFinder::addNewSearchResult(const std::list<pair<int, int> > &newlines) 688 void SDLineFinder::addNewSearchResult(const std::list<pair<int, int> > &newlines, 689 std::list<std::pair<int, int> > &lines_list) 572 690 throw(AipsError) 573 691 { … … 578 696 // the first item, which has a non-void intersection or touches 579 697 // the new line 580 std::list<pair<int,int> >::iterator pos_beg=find_if(lines .begin(),581 lines .end(), IntersectsWith(*cli));698 std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(), 699 lines_list.end(), IntersectsWith(*cli)); 582 700 // the last such item 583 701 std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg, 584 lines .end(), not1(IntersectsWith(*cli)));702 lines_list.end(), not1(IntersectsWith(*cli))); 585 703 586 704 // extract all lines which intersect or touch a new one into … … 589 707 // one. 590 708 std::list<pair<int,int> > lines_buffer; 591 lines_buffer.splice(lines_buffer.end(),lines , pos_beg, pos_end);709 lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end); 592 710 593 711 // build a union of all intersecting lines … … 596 714 597 715 // search for a right place for the new line (union_line) and add 598 std::list<pair<int,int> >::iterator pos2insert=find_if(lines .begin(),599 lines .end(), LaterThan(union_line));600 lines .insert(pos2insert,union_line);716 std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(), 717 lines_list.end(), LaterThan(union_line)); 718 lines_list.insert(pos2insert,union_line); 601 719 } 602 720 } … … 608 726 } 609 727 } 728 729 // extend all line ranges to the point where a value stored in the 730 // specified vector changes (e.g. value-mean change its sign) 731 // This operation is necessary to include line wings, which are below 732 // the detection threshold. If lines becomes adjacent, they are 733 // merged together. Any masked channel stops the extension 734 void SDLineFinder::searchForWings(std::list<std::pair<int, int> > &newlines, 735 const casa::Vector<casa::Int> &signs) throw(casa::AipsError) 736 { 737 try { 738 for (std::list<pair<int,int> >::iterator li=newlines.begin(); 739 li!=newlines.end();++li) { 740 // update the left hand side 741 for (int n=li->first-1;n>=edge.first;--n) { 742 if (!mask[n]) break; 743 if (signs[n]==signs[li->first] && signs[li->first]) 744 li->first=n; 745 else break; 746 } 747 // update the right hand side 748 for (int n=li->second;n<edge.second;++n) { 749 if (!mask[n]) break; 750 if (signs[n]==signs[li->second-1] && signs[li->second-1]) 751 li->second=n; 752 else break; 753 } 754 } 755 // need to search for possible mergers. 756 std::list<std::pair<int, int> > result_buffer; 757 addNewSearchResult(newlines,result_buffer); 758 newlines.clear(); 759 newlines.splice(newlines.end(),result_buffer); 760 } 761 catch (const AipsError &ae) { 762 throw; 763 } 764 catch (const exception &ex) { 765 throw AipsError(String("SDLineFinder::extendLines - STL error: ")+ex.what()); 766 } 767 } -
trunk/src/SDLineFinder.h
r343 r344 89 89 // be adjacent or have a non-void intersection, they are joined into 90 90 // the new line 91 void addNewSearchResult(const std::list<std::pair<int, int> > &newlines) 91 static void addNewSearchResult(const std::list<std::pair<int, int> > 92 &newlines, std::list<std::pair<int, int> > &lines_list) 92 93 throw(casa::AipsError); 94 95 // extend all line ranges to the point where a value stored in the 96 // specified vector changes (e.g. value-mean change its sign) 97 // This operation is necessary to include line wings, which are below 98 // the detection threshold. If lines becomes adjacent, they are 99 // merged together. Any masked channel stops the extension 100 void searchForWings(std::list<std::pair<int, int> > &newlines, 101 const casa::Vector<casa::Int> &signs) 102 throw(casa::AipsError); 93 103 94 104 // An auxiliary object function to test whether two lines have a non-void
Note:
See TracChangeset
for help on using the changeset viewer.