Changeset 996 for trunk/src/STLineFinder.cpp
- Timestamp:
- 04/06/06 13:45:58 (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/STLineFinder.cpp
r991 r996 81 81 int start_advance; // number of channel from which the box can 82 82 // be moved (the middle of the box, if there is no 83 83 // masking) 84 84 public: 85 85 // set up the object with the references to actual data … … 87 87 RunningBox(const casa::Vector<casa::Float> &in_spectrum, 88 88 const casa::Vector<casa::Bool> &in_mask, 89 90 89 const std::pair<int,int> &in_edge, 90 int in_max_box_nchan) throw(AipsError); 91 91 92 92 // access to the statistics … … 139 139 int min_nchan; // A minimum number of consequtive 140 140 // channels, which should satisfy 141 142 141 // the detection criterion, to be 142 // a detection 143 143 casa::Float threshold; // detection threshold - the 144 144 // minimal signal to noise ratio … … 148 148 casa::Vector<Int> signs; // An array to store the signs of 149 149 // the value - current mean 150 150 // (used to search wings) 151 151 casa::Int last_sign; // a sign (+1, -1 or 0) of the 152 152 // last point of the detected line … … 157 157 LFAboveThreshold(std::list<pair<int,int> > &in_lines, 158 158 int in_min_nchan = 3, 159 159 casa::Float in_threshold = 5) throw(); 160 160 virtual ~LFAboveThreshold() throw(); 161 161 … … 175 175 // max_box_nchan - number of channels in the running box 176 176 void findLines(const casa::Vector<casa::Float> &spectrum, 177 178 179 177 const casa::Vector<casa::Bool> &mask, 178 const std::pair<int,int> &edge, 179 int max_box_nchan) throw(casa::AipsError); 180 180 181 181 protected: … … 213 213 RunningBox::RunningBox(const casa::Vector<casa::Float> &in_spectrum, 214 214 const casa::Vector<casa::Bool> &in_mask, 215 216 215 const std::pair<int,int> &in_edge, 216 int in_max_box_nchan) throw(AipsError) : 217 217 spectrum(in_spectrum), mask(in_mask), edge(in_edge), 218 218 max_box_nchan(in_max_box_nchan) 219 219 { 220 220 rewind(); … … 287 287 sumf+=spectrum[ch]; 288 288 sumf2+=square(spectrum[ch]); 289 290 291 292 289 sumch+=Float(ch); 290 sumch2+=square(Float(ch)); 291 sumfch+=spectrum[ch]*Float(ch); 292 need2recalculate=True; 293 293 } 294 294 int ch2remove=ch-max_box_nchan; … … 298 298 sumf-=spectrum[ch2remove]; 299 299 sumf2-=square(spectrum[ch2remove]); 300 301 302 303 300 sumch-=Float(ch2remove); 301 sumch2-=square(Float(ch2remove)); 302 sumfch-=spectrum[ch2remove]*Float(ch2remove); 303 need2recalculate=True; 304 304 } 305 305 } … … 361 361 casa::Float in_threshold) throw() : 362 362 min_nchan(in_min_nchan), threshold(in_threshold), 363 363 lines(in_lines), running_box(NULL) {} 364 364 365 365 LFAboveThreshold::~LFAboveThreshold() throw() … … 427 427 if (is_detected_before) { 428 428 if (cur_line.second-cur_line.first>min_nchan) { 429 430 431 432 433 434 435 436 437 438 439 440 429 // it was a detection. We need to change the list 430 Bool add_new_line=False; 431 if (lines.size()) { 432 for (int i=lines.back().second;i<cur_line.first;++i) 433 if (mask[i]) { // one valid channel in between 434 // means that we deal with a separate line 435 add_new_line=True; 436 break; 437 } 438 } else add_new_line=True; 439 if (add_new_line) 440 lines.push_back(cur_line); 441 441 else lines.back().second=cur_line.second; 442 443 442 } 443 is_detected_before=False; 444 444 } 445 445 } … … 463 463 // find spectral lines and add them into list 464 464 void LFAboveThreshold::findLines(const casa::Vector<casa::Float> &spectrum, 465 466 467 465 const casa::Vector<casa::Bool> &mask, 466 const std::pair<int,int> &edge, 467 int max_box_nchan) 468 468 throw(casa::AipsError) 469 469 { … … 483 483 for (;running_box->haveMore();running_box->next()) 484 484 variances[running_box->getChannel()-edge.first]= 485 485 running_box->getLinVariance(); 486 486 487 487 // in the future we probably should do a proper Chi^2 estimation … … 511 511 const int ch=running_box->getChannel(); 512 512 if (running_box->getNumberOfBoxPoints()>=minboxnchan) 513 514 515 516 517 518 513 processChannel(mask[ch] && (fabs(running_box->aboveMean()) >= 514 threshold*offline_variance), mask); 515 else processCurLine(mask); // just finish what was accumulated before 516 517 signs[ch]=getAboveMeanSign(); 518 // os<<ch<<" "<<spectrum[ch]<<" "<<fabs(running_box->aboveMean())<<" "<< 519 519 // threshold*offline_variance<<endl; 520 520 521 522 523 524 525 526 521 const Float buf=running_box->aboveMean(); 522 if (buf>0) signs[ch]=1; 523 else if (buf<0) signs[ch]=-1; 524 else if (buf==0) signs[ch]=0; 525 // os<<ch<<" "<<spectrum[ch]<<" "<<running_box->getLinMean()<<" "<< 526 // threshold*offline_variance<<endl; 527 527 } 528 528 if (lines.size()) … … 653 653 void STLineFinder::setOptions(const casa::Float &in_threshold, 654 654 const casa::Int &in_min_nchan, 655 655 const casa::Int &in_avg_limit, 656 656 const casa::Float &in_box_size) throw() 657 657 { … … 680 680 // can be achieved using a spectrum mask only 681 681 int STLineFinder::findLines(const std::vector<bool> &in_mask, 682 683 684 { 685 const int minboxnchan=4;682 const std::vector<int> &in_edge, 683 const casa::uInt &whichRow) throw(casa::AipsError) 684 { 685 //const int minboxnchan=4; 686 686 if (scan.null()) 687 687 throw AipsError("STLineFinder::findLines - a scan should be set first," … … 704 704 if (in_edge.size()>2) 705 705 throw AipsError("STLineFinder::findLines - the length of the in_edge parameter" 706 706 "should not exceed 2"); 707 707 if (!in_edge.size()) { 708 708 // all spectra, no rejection 709 709 edge.first=0; 710 710 edge.second=nchan; 711 711 } else { 712 712 edge.first=in_edge[0]; 713 714 715 716 if (edge.first>= nchan)717 713 if (edge.first<0) 714 throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative" 715 "number of channels to drop"); 716 if (edge.first>=int(nchan)) 717 throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter"); 718 718 if (in_edge.size()==2) { 719 720 721 722 719 edge.second=in_edge[1]; 720 if (edge.second<0) 721 throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative" 722 "number of channels to drop"); 723 723 edge.second=nchan-edge.second; 724 724 } else edge.second=nchan-edge.first; 725 725 if (edge.second<0 || (edge.first>=edge.second)) 726 726 throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter"); 727 727 } 728 728 … … 743 743 Int avg_factor=1; // this number of adjacent channels is averaged together 744 744 // the total number of the channels is not altered 745 746 745 // instead, min_nchan is also scaled 746 // it helps to search for broad lines 747 747 Vector<Int> signs; // a buffer for signs of the value - mean quantity 748 748 // see LFAboveThreshold for details 749 750 751 749 // We need only signs resulted from last iteration 750 // because all previous values may be corrupted by the 751 // presence of spectral lines 752 752 while (true) { 753 753 // a buffer for new lines found at this iteration … … 758 758 LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold); 759 759 lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan); 760 761 760 signs.resize(lfalg.getSigns().nelements()); 761 signs=lfalg.getSigns(); 762 762 first_pass=False; 763 763 if (!new_lines.size()) 764 765 764 throw AipsError("spurious"); // nothing new - use the same 765 // code as for a real exception 766 766 } 767 767 catch(const AipsError &ae) { 768 768 if (first_pass) throw; 769 769 // nothing new - proceed to the next step of averaging, if any 770 771 772 773 774 775 776 777 770 // (to search for broad lines) 771 if (avg_factor>avg_limit) break; // averaging up to avg_limit 772 // adjacent channels, 773 // stop after that 774 avg_factor*=2; // twice as more averaging 775 subtractBaseline(temp_mask,9); 776 averageAdjacentChannels(temp_mask,avg_factor); 777 continue; 778 778 } 779 779 keepStrongestOnly(temp_mask,new_lines,max_box_nchan); … … 804 804 Fitter sdf; 805 805 std::vector<float> absc(spectrum.nelements()); 806 for ( Int i=0;i<absc.size();++i)806 for (unsigned int i=0;i<absc.size();++i) 807 807 absc[i]=float(i)/float(spectrum.nelements()); 808 808 std::vector<float> spec; … … 835 835 for (int k=n;k<n+boxsize && k<edge.second;++k) 836 836 if (mask2update[k]) { // k is a valid channel 837 838 837 mean+=spectrum[k]; 838 ++nboxch; 839 839 } 840 840 if (nboxch<boxsize) // mask these channels 841 841 for (int k=n;k<n+boxsize && k<edge.second;++k) 842 842 mask2update[k]=False; 843 843 else { 844 844 mean/=Float(boxsize); 845 846 845 for (int k=n;k<n+boxsize && k<edge.second;++k) 846 spectrum[k]=mean; 847 847 } 848 848 } … … 866 866 if (!lines.size()) 867 867 throw AipsError("STLineFinder::getMask - one have to search for " 868 868 "lines first, use find_lines"); 869 869 */ 870 870 std::vector<bool> res_mask(mask.nelements()); 871 871 // iterator through lines 872 872 std::list<std::pair<int,int> >::const_iterator cli=lines.begin(); 873 for (int ch=0;ch< res_mask.size();++ch)873 for (int ch=0;ch<int(res_mask.size());++ch) 874 874 if (ch<edge.first || ch>=edge.second) res_mask[ch]=false; 875 876 877 878 879 880 875 else if (!mask[ch]) res_mask[ch]=false; 876 else { 877 res_mask[ch]=!invert; // no line by default 878 if (cli==lines.end()) continue; 879 if (ch>=cli->first && ch<cli->second) 880 res_mask[ch]=invert; // this is a line 881 881 if (ch>=cli->second) 882 883 882 ++cli; // next line in the list 883 } 884 884 885 885 return res_mask; … … 925 925 if (!lines.size()) 926 926 throw AipsError("STLineFinder::getLineRangesInChannels - one have to search for " 927 927 "lines first, use find_lines"); 928 928 929 929 std::vector<int> res(2*lines.size()); … … 932 932 std::vector<int>::iterator ri=res.begin(); 933 933 for (;cli!=lines.end() && ri!=res.end();++cli,++ri) { 934 935 936 934 *ri=cli->first; 935 if (++ri!=res.end()) 936 *ri=cli->second-1; 937 937 } 938 938 return res; … … 962 962 // max_box_nchan - channels in the running box for baseline filtering 963 963 void STLineFinder::keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask, 964 965 964 std::list<std::pair<int, int> > &lines2update, 965 int max_box_nchan) 966 966 throw (casa::AipsError) 967 967 { … … 982 982 983 983 if (li==lines2update.end()) break; // no more lines 984 985 986 987 988 989 990 991 992 993 994 995 984 const int ch=running_box.getChannel(); 985 if (ch>=li->first && ch<li->second) 986 if (temp_mask[ch] && tmp_flux<fabs(running_box.aboveMean())) 987 tmp_flux=fabs(running_box.aboveMean()); 988 if (ch==li->second-1) { 989 if (peak_flux<tmp_flux) { // if peak_flux=-1, this condition 990 peak_flux=tmp_flux; // will be satisfied 991 strongli=li; 992 } 993 ++li; 994 tmp_flux=-1; 995 } 996 996 } 997 997 std::list<std::pair<int,int> > res; … … 1029 1029 cli!=newlines.end();++cli) { 1030 1030 1031 1032 1033 1034 1035 1036 1037 1031 // the first item, which has a non-void intersection or touches 1032 // the new line 1033 std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(), 1034 lines_list.end(), IntersectsWith(*cli)); 1035 // the last such item 1036 std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg, 1037 lines_list.end(), not1(IntersectsWith(*cli))); 1038 1038 1039 1039 // extract all lines which intersect or touch a new one into 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1040 // a temporary buffer. This may invalidate the iterators 1041 // line_buffer may be empty, if no lines intersects with a new 1042 // one. 1043 std::list<pair<int,int> > lines_buffer; 1044 lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end); 1045 1046 // build a union of all intersecting lines 1047 pair<int,int> union_line=for_each(lines_buffer.begin(), 1048 lines_buffer.end(),BuildUnion(*cli)).result(); 1049 1050 // search for a right place for the new line (union_line) and add 1051 std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(), 1052 lines_list.end(), LaterThan(union_line)); 1053 lines_list.insert(pos2insert,union_line); 1054 1054 } 1055 1055 } … … 1069 1069 void LFLineListOperations::searchForWings(std::list<std::pair<int, int> > &newlines, 1070 1070 const casa::Vector<casa::Int> &signs, 1071 1072 1071 const casa::Vector<casa::Bool> &mask, 1072 const std::pair<int,int> &edge) throw(casa::AipsError) 1073 1073 { 1074 1074 try { 1075 1075 for (std::list<pair<int,int> >::iterator li=newlines.begin(); 1076 1076 li!=newlines.end();++li) { 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1077 // update the left hand side 1078 for (int n=li->first-1;n>=edge.first;--n) { 1079 if (!mask[n]) break; 1080 if (signs[n]==signs[li->first] && signs[li->first]) 1081 li->first=n; 1082 else break; 1083 } 1084 // update the right hand side 1085 for (int n=li->second;n<edge.second;++n) { 1086 if (!mask[n]) break; 1087 if (signs[n]==signs[li->second-1] && signs[li->second-1]) 1088 li->second=n; 1089 else break; 1090 } 1091 1091 } 1092 1092 // need to search for possible mergers.
Note: See TracChangeset
for help on using the changeset viewer.