- Timestamp:
- 03/08/06 11:58:43 (19 years ago)
- Location:
- trunk/src
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/SDLineFinder.cc
r551 r881 1 1 //#--------------------------------------------------------------------------- 2 //# S DLineFinder.cc: A class for automated spectral line search2 //# STLineFinder.cc: A class for automated spectral line search 3 3 //#-------------------------------------------------------------------------- 4 4 //# Copyright (C) 2004 … … 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id: 29 //# $Id:$ 30 30 //#--------------------------------------------------------------------------- 31 31 … … 50 50 /////////////////////////////////////////////////////////////////////////////// 51 51 // 52 // RunningBox - a running box calculator. This class implements 52 // RunningBox - a running box calculator. This class implements 53 53 // interations over the specified spectrum and calculates 54 54 // running box filter statistics. … … 57 57 class RunningBox { 58 58 // The input data to work with. Use reference symantics to avoid 59 // an unnecessary copying 59 // an unnecessary copying 60 60 const casa::Vector<casa::Float> &spectrum; // a buffer for the spectrum 61 61 const casa::Vector<casa::Bool> &mask; // associated mask 62 62 const std::pair<int,int> &edge; // start and stop+1 channels 63 63 // to work with 64 64 65 65 // statistics for running box filtering 66 66 casa::Float sumf; // sum of fluxes … … 69 69 casa::Float sumch2; // sum of squares of channel numbers (for linear fit) 70 70 casa::Float sumfch; // sum of flux*(channel number) (for linear fit) 71 71 72 72 int box_chan_cntr; // actual number of channels in the box 73 73 int max_box_nchan; // maximum allowed number of channels in the box … … 90 90 const std::pair<int,int> &in_edge, 91 91 int in_max_box_nchan) throw(AipsError); 92 92 93 93 // access to the statistics 94 94 const casa::Float& getLinMean() const throw(AipsError); … … 96 96 const casa::Float aboveMean() const throw(AipsError); 97 97 int getChannel() const throw(); 98 98 99 99 // actual number of channels in the box (max_box_nchan, if no channels 100 100 // are masked) … … 109 109 // go to start 110 110 void rewind() throw(AipsError); 111 111 112 112 protected: 113 113 // supplementary function to control running mean calculations. … … 142 142 // the detection criterion, to be 143 143 // a detection 144 casa::Float threshold; // detection threshold - the 144 casa::Float threshold; // detection threshold - the 145 145 // minimal signal to noise ratio 146 146 std::list<pair<int,int> > &lines; // list where detections are saved … … 157 157 casa::Float in_threshold = 5) throw(); 158 158 virtual ~LFAboveThreshold() throw(); 159 159 160 160 // replace the detection criterion 161 161 void setCriterion(int in_min_nchan, casa::Float in_threshold) throw(); … … 198 198 /////////////////////////////////////////////////////////////////////////////// 199 199 // 200 // RunningBox - a running box calculator. This class implements 200 // RunningBox - a running box calculator. This class implements 201 201 // interations over the specified spectrum and calculates 202 202 // running box filter statistics. … … 227 227 ++initial_box_ch) 228 228 advanceRunningBox(initial_box_ch); 229 230 if (initial_box_ch==edge.second) 229 230 if (initial_box_ch==edge.second) 231 231 throw AipsError("RunningBox::rewind - too much channels are masked"); 232 232 233 233 cur_channel=edge.first; 234 start_advance=initial_box_ch-max_box_nchan/2; 234 start_advance=initial_box_ch-max_box_nchan/2; 235 235 } 236 236 … … 291 291 --box_chan_cntr; 292 292 sumf-=spectrum[ch2remove]; 293 sumf2-=square(spectrum[ch2remove]); 293 sumf2-=square(spectrum[ch2remove]); 294 294 sumch-=Float(ch2remove); 295 295 sumch2-=square(Float(ch2remove)); … … 319 319 { 320 320 AlwaysAssert(box_chan_cntr, AipsError); 321 321 322 322 Float mean=sumf/Float(box_chan_cntr); 323 323 … … 385 385 cur_line.second=running_box->getChannel()+1; 386 386 } 387 } else processCurLine(mask); 387 } else processCurLine(mask); 388 388 } 389 389 catch (const AipsError &ae) { 390 390 throw; 391 } 391 } 392 392 catch (const exception &ex) { 393 393 throw AipsError(String("LFAboveThreshold::processChannel - STL error: ")+ex.what()); … … 401 401 { 402 402 try { 403 if (is_detected_before) { 403 if (is_detected_before) { 404 404 if (cur_line.second-cur_line.first>min_nchan) { 405 405 // it was a detection. We need to change the list 406 406 Bool add_new_line=False; 407 if (lines.size()) { 407 if (lines.size()) { 408 408 for (int i=lines.back().second;i<cur_line.first;++i) 409 409 if (mask[i]) { // one valid channel in between … … 412 412 break; 413 413 } 414 } else add_new_line=True; 415 if (add_new_line) 414 } else add_new_line=True; 415 if (add_new_line) 416 416 lines.push_back(cur_line); 417 else lines.back().second=cur_line.second; 417 else lines.back().second=cur_line.second; 418 418 } 419 419 is_detected_before=False; 420 } 420 } 421 421 } 422 422 catch (const AipsError &ae) { 423 423 throw; 424 } 424 } 425 425 catch (const exception &ex) { 426 426 throw AipsError(String("LFAboveThreshold::processCurLine - STL error: ")+ex.what()); … … 450 450 running_box=new RunningBox(spectrum,mask,edge,max_box_nchan); 451 451 452 452 453 453 // determine the off-line variance first 454 454 // an assumption made: lines occupy a small part of the spectrum 455 455 456 456 std::vector<float> variances(edge.second-edge.first); 457 457 DebugAssert(variances.size(),AipsError); 458 459 for (;running_box->haveMore();running_box->next()) 458 459 for (;running_box->haveMore();running_box->next()) 460 460 variances[running_box->getChannel()-edge.first]= 461 461 running_box->getLinVariance(); 462 462 463 463 // in the future we probably should do a proper Chi^2 estimation 464 464 // now a simple 80% of smaller values will be used. … … 466 466 // due to a bias of the Chi^2 distribution. 467 467 stable_sort(variances.begin(),variances.end()); 468 468 469 469 Float offline_variance=0; 470 uInt offline_cnt=uInt(0.8*variances.size()); 470 uInt offline_cnt=uInt(0.8*variances.size()); 471 471 if (!offline_cnt) offline_cnt=variances.size(); // no much else left, 472 472 // although it is very inaccurate 473 for (uInt n=0;n<offline_cnt;++n) 474 offline_variance+=variances[n]; 475 offline_variance/=Float(offline_cnt); 476 473 for (uInt n=0;n<offline_cnt;++n) 474 offline_variance+=variances[n]; 475 offline_variance/=Float(offline_cnt); 476 477 477 // actual search algorithm 478 478 is_detected_before=False; … … 502 502 catch (const AipsError &ae) { 503 503 throw; 504 } 504 } 505 505 catch (const exception &ex) { 506 506 throw AipsError(String("LFAboveThreshold::findLines - STL error: ")+ex.what()); … … 583 583 if (line2.second<line1.first) return false; // line2 is at lower channels 584 584 if (line2.first>line1.second) return true; // line2 is at upper channels 585 585 586 586 // line2 intersects with line1. We should have no such situation in 587 587 // practice … … 595 595 /////////////////////////////////////////////////////////////////////////////// 596 596 // 597 // S DLineFinder - a class for automated spectral line search598 // 599 // 600 601 S DLineFinder::SDLineFinder() throw() : edge(0,0)597 // STLineFinder - a class for automated spectral line search 598 // 599 // 600 601 STLineFinder::STLineFinder() throw() : edge(0,0) 602 602 { 603 603 setOptions(); … … 614 614 // in_avg_limit perform the averaging of no more than in_avg_limit 615 615 // adjacent channels to search for broad lines 616 // Default is 8, but for a bad baseline shape this 616 // Default is 8, but for a bad baseline shape this 617 617 // parameter should be decreased (may be even down to a 618 618 // minimum of 1 to disable this option) to avoid 619 619 // confusing of baseline undulations with a real line. 620 // Setting a very large value doesn't usually provide 621 // valid detections. 620 // Setting a very large value doesn't usually provide 621 // valid detections. 622 622 // in_box_size the box size for running mean calculation. Default is 623 623 // 1./5. of the whole spectrum size 624 void S DLineFinder::setOptions(const casa::Float &in_threshold,624 void STLineFinder::setOptions(const casa::Float &in_threshold, 625 625 const casa::Int &in_min_nchan, 626 626 const casa::Int &in_avg_limit, … … 633 633 } 634 634 635 S DLineFinder::~SDLineFinder() throw(AipsError) {}635 STLineFinder::~STLineFinder() throw(AipsError) {} 636 636 637 637 // set scan to work with (in_scan parameter), associated mask (in_mask … … 641 641 // channels to drop from both sides of the spectrum 642 642 // in_edge is introduced for convinience, although all functionality 643 // can be achieved using a spectrum mask only 644 void S DLineFinder::setScan(const SDMemTableWrapper &in_scan,643 // can be achieved using a spectrum mask only 644 void STLineFinder::setScan(const ScantableWrapper &in_scan, 645 645 const std::vector<bool> &in_mask, 646 646 const boost::python::tuple &in_edge) throw(AipsError) … … 651 651 652 652 mask=in_mask; 653 if (mask.nelements()!=scan->n Chan())654 throw AipsError("S DLineFinder::setScan - in_scan and in_mask have different"653 if (mask.nelements()!=scan->nchan()) 654 throw AipsError("STLineFinder::setScan - in_scan and in_mask have different" 655 655 "number of spectral channels."); 656 656 … … 658 658 int n=extract<int>(in_edge.attr("__len__")()); 659 659 if (n>2 || n<0) 660 throw AipsError("S DLineFinder::setScan - the length of the in_edge parameter"660 throw AipsError("STLineFinder::setScan - the length of the in_edge parameter" 661 661 "should not exceed 2"); 662 662 if (!n) { 663 // all spectr um, no rejection663 // all spectra, no rejection 664 664 edge.first=0; 665 edge.second=scan->n Chan();665 edge.second=scan->nchan(); 666 666 } else { 667 667 edge.first=extract<int>(in_edge.attr("__getitem__")(0)); 668 668 if (edge.first<0) 669 throw AipsError("S DLineFinder::setScan - the in_edge parameter has a negative"669 throw AipsError("STLineFinder::setScan - the in_edge parameter has a negative" 670 670 "number of channels to drop"); 671 if (edge.first>=scan->n Chan())672 throw AipsError("S DLineFinder::setScan - all channels are rejected by the in_edge parameter");671 if (edge.first>=scan->nchan()) 672 throw AipsError("STLineFinder::setScan - all channels are rejected by the in_edge parameter"); 673 673 if (n==2) { 674 674 edge.second=extract<int>(in_edge.attr("__getitem__")(1)); 675 675 if (edge.second<0) 676 throw AipsError("S DLineFinder::setScan - the in_edge parameter has a negative"676 throw AipsError("STLineFinder::setScan - the in_edge parameter has a negative" 677 677 "number of channels to drop"); 678 edge.second=scan->n Chan()-edge.second;679 } else edge.second=scan->n Chan()-edge.first;678 edge.second=scan->nchan()-edge.second; 679 } else edge.second=scan->nchan()-edge.first; 680 680 if (edge.second<0 || (edge.first>=edge.second)) 681 throw AipsError("S DLineFinder::setScan - all channels are rejected by the in_edge parameter");682 } 681 throw AipsError("STLineFinder::setScan - all channels are rejected by the in_edge parameter"); 682 } 683 683 } 684 684 catch(const AipsError &ae) { 685 685 // setScan is unsuccessfull, reset scan/mask/edge 686 scan=CountedConstPtr<S DMemTable>(); // null pointer686 scan=CountedConstPtr<Scantable>(); // null pointer 687 687 mask.resize(0); 688 688 edge=pair<int,int>(0,0); … … 692 692 693 693 // search for spectral lines. Number of lines found is returned 694 int S DLineFinder::findLines(const casa::uInt &whichRow) throw(casa::AipsError)695 { 694 int STLineFinder::findLines(const casa::uInt &whichRow) throw(casa::AipsError) 695 { 696 696 const int minboxnchan=4; 697 697 if (scan.null()) 698 throw AipsError("S DLineFinder::findLines - a scan should be set first,"698 throw AipsError("STLineFinder::findLines - a scan should be set first," 699 699 " use set_scan"); 700 DebugAssert(mask.nelements()==scan->n Chan(), AipsError);701 int max_box_nchan=int(scan->n Chan()*box_size); // number of channels in running700 DebugAssert(mask.nelements()==scan->nchan(), AipsError); 701 int max_box_nchan=int(scan->nchan()*box_size); // number of channels in running 702 702 // box 703 703 if (max_box_nchan<2) 704 throw AipsError("SDLineFinder::findLines - box_size is too small"); 705 706 scan->getSpectrum(spectrum, whichRow); 704 throw AipsError("STLineFinder::findLines - box_size is too small"); 705 706 spectrum.resize(); 707 spectrum = Vector<Float>(scan->getSpectrum(whichRow)); 707 708 708 709 lines.resize(0); // search from the scratch … … 722 723 while (true) { 723 724 // a buffer for new lines found at this iteration 724 std::list<pair<int,int> > new_lines; 725 std::list<pair<int,int> > new_lines; 725 726 726 727 try { … … 745 746 subtractBaseline(temp_mask,9); 746 747 averageAdjacentChannels(temp_mask,avg_factor); 747 continue; 748 continue; 748 749 } 749 750 keepStrongestOnly(temp_mask,new_lines,max_box_nchan); … … 751 752 addNewSearchResult(new_lines,lines); 752 753 // get a new mask 753 temp_mask=getMask(); 754 } 755 754 temp_mask=getMask(); 755 } 756 756 757 // an additional search for wings because in the presence of very strong 757 758 // lines temporary mean used at each iteration will be higher than 758 759 // the true mean 759 760 760 761 if (lines.size()) 761 762 LFLineListOperations::searchForWings(lines,signs,mask,edge); 762 763 763 764 return int(lines.size()); 764 765 } … … 767 768 // spectrum. It uses the SDFitter class. This action is required before 768 769 // reducing the spectral resolution if the baseline shape is bad 769 void S DLineFinder::subtractBaseline(const casa::Vector<casa::Bool> &temp_mask,770 void STLineFinder::subtractBaseline(const casa::Vector<casa::Bool> &temp_mask, 770 771 const casa::Int &order) throw(casa::AipsError) 771 772 { … … 783 784 sdf.setExpression("poly",order); 784 785 if (!sdf.fit()) return; // fit failed, use old spectrum 785 spectrum=casa::Vector<casa::Float>(sdf.getResidual()); 786 spectrum=casa::Vector<casa::Float>(sdf.getResidual()); 786 787 } 787 788 … … 792 793 // array specified, rather than the field of this class 793 794 // boxsize - a number of adjacent channels to average 794 void S DLineFinder::averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update,795 void STLineFinder::averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update, 795 796 const casa::Int &boxsize) 796 797 throw(casa::AipsError) … … 798 799 DebugAssert(mask2update.nelements()==spectrum.nelements(), AipsError); 799 800 DebugAssert(boxsize!=0,AipsError); 800 801 801 802 for (int n=edge.first;n<edge.second;n+=boxsize) { 802 803 DebugAssert(n<spectrum.nelements(),AipsError); … … 807 808 mean+=spectrum[k]; 808 809 ++nboxch; 809 } 810 } 810 811 if (nboxch<boxsize) // mask these channels 811 812 for (int k=n;k<n+boxsize && k<edge.second;++k) … … 825 826 // in setScan) or dropped out by the edge parameter (in_edge 826 827 // in setScan) are still excluded regardless on the invert option 827 std::vector<bool> S DLineFinder::getMask(bool invert)828 std::vector<bool> STLineFinder::getMask(bool invert) 828 829 const throw(casa::AipsError) 829 830 { 830 831 try { 831 832 if (scan.null()) 832 throw AipsError("S DLineFinder::getMask - a scan should be set first,"833 throw AipsError("STLineFinder::getMask - a scan should be set first," 833 834 " use set_scan followed by find_lines"); 834 DebugAssert(mask.nelements()==scan->n Chan(), AipsError);835 DebugAssert(mask.nelements()==scan->nchan(), AipsError); 835 836 /* 836 837 if (!lines.size()) 837 throw AipsError("S DLineFinder::getMask - one have to search for "838 throw AipsError("STLineFinder::getMask - one have to search for " 838 839 "lines first, use find_lines"); 839 */ 840 */ 840 841 std::vector<bool> res_mask(mask.nelements()); 841 842 // iterator through lines 842 843 std::list<std::pair<int,int> >::const_iterator cli=lines.begin(); 843 for (int ch=0;ch<res_mask.size();++ch) 844 for (int ch=0;ch<res_mask.size();++ch) 844 845 if (ch<edge.first || ch>=edge.second) res_mask[ch]=false; 845 846 else if (!mask[ch]) res_mask[ch]=false; 846 else { 847 else { 847 848 res_mask[ch]=!invert; // no line by default 848 849 if (cli==lines.end()) continue; … … 852 853 ++cli; // next line in the list 853 854 } 854 855 855 856 return res_mask; 856 857 } 857 858 catch (const AipsError &ae) { 858 859 throw; 859 } 860 } 860 861 catch (const exception &ex) { 861 throw AipsError(String("S DLineFinder::getMask - STL error: ")+ex.what());862 throw AipsError(String("STLineFinder::getMask - STL error: ")+ex.what()); 862 863 } 863 864 } … … 865 866 // get range for all lines found. The same units as used in the scan 866 867 // will be returned (e.g. velocity instead of channels). 867 std::vector<double> S DLineFinder::getLineRanges()868 std::vector<double> STLineFinder::getLineRanges() 868 869 const throw(casa::AipsError) 869 870 { … … 877 878 for (;cri!=ranges.end() && outi!=res.end();++cri,++outi) 878 879 if (uInt(*cri)>=vel.size()) 879 throw AipsError("S DLineFinder::getLineRanges - getAbcissa provided less channels than reqired");880 throw AipsError("STLineFinder::getLineRanges - getAbcissa provided less channels than reqired"); 880 881 else *outi=vel[*cri]; 881 882 return res; … … 884 885 // The same as getLineRanges, but channels are always used to specify 885 886 // the range 886 std::vector<int> S DLineFinder::getLineRangesInChannels()887 std::vector<int> STLineFinder::getLineRangesInChannels() 887 888 const throw(casa::AipsError) 888 889 { 889 890 try { 890 891 if (scan.null()) 891 throw AipsError("S DLineFinder::getLineRangesInChannels - a scan should be set first,"892 throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first," 892 893 " use set_scan followed by find_lines"); 893 DebugAssert(mask.nelements()==scan->n Chan(), AipsError);894 894 DebugAssert(mask.nelements()==scan->nchan(), AipsError); 895 895 896 if (!lines.size()) 896 throw AipsError("S DLineFinder::getLineRangesInChannels - one have to search for "897 throw AipsError("STLineFinder::getLineRangesInChannels - one have to search for " 897 898 "lines first, use find_lines"); 898 899 899 900 std::vector<int> res(2*lines.size()); 900 901 // iterator through lines & result 901 902 std::list<std::pair<int,int> >::const_iterator cli=lines.begin(); 902 903 std::vector<int>::iterator ri=res.begin(); 903 for (;cli!=lines.end() && ri!=res.end();++cli,++ri) { 904 for (;cli!=lines.end() && ri!=res.end();++cli,++ri) { 904 905 *ri=cli->first; 905 if (++ri!=res.end()) 906 *ri=cli->second-1; 907 } 906 if (++ri!=res.end()) 907 *ri=cli->second-1; 908 } 908 909 return res; 909 910 } 910 911 catch (const AipsError &ae) { 911 912 throw; 912 } 913 } 913 914 catch (const exception &ex) { 914 throw AipsError(String("S DLineFinder::getLineRanges - STL error: ")+ex.what());915 throw AipsError(String("STLineFinder::getLineRanges - STL error: ")+ex.what()); 915 916 } 916 917 } … … 920 921 // an auxiliary function to remove all lines from the list, except the 921 922 // strongest one (by absolute value). If the lines removed are real, 922 // they will be find again at the next iteration. This approach 923 // increases the number of iterations required, but is able to remove 923 // they will be find again at the next iteration. This approach 924 // increases the number of iterations required, but is able to remove 924 925 // the sidelobes likely to occur near strong lines. 925 926 // Later a better criterion may be implemented, e.g. 926 927 // taking into consideration the brightness of different lines. Now 927 // use the simplest solution 928 // use the simplest solution 928 929 // temp_mask - mask to work with (may be different from original mask as 929 930 // the lines previously found may be masked) … … 931 932 // nothing will be done if it is empty 932 933 // max_box_nchan - channels in the running box for baseline filtering 933 void S DLineFinder::keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask,934 void STLineFinder::keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask, 934 935 std::list<std::pair<int, int> > &lines2update, 935 936 int max_box_nchan) … … 952 953 953 954 if (li==lines2update.end()) break; // no more lines 954 const int ch=running_box.getChannel(); 955 const int ch=running_box.getChannel(); 955 956 if (ch>=li->first && ch<li->second) 956 957 if (temp_mask[ch] && tmp_flux<fabs(running_box.aboveMean())) … … 964 965 tmp_flux=-1; 965 966 } 966 } 967 } 967 968 std::list<std::pair<int,int> > res; 968 969 res.splice(res.end(),lines2update,strongli); … … 972 973 catch (const AipsError &ae) { 973 974 throw; 974 } 975 } 975 976 catch (const exception &ex) { 976 throw AipsError(String("S DLineFinder::keepStrongestOnly - STL error: ")+ex.what());977 throw AipsError(String("STLineFinder::keepStrongestOnly - STL error: ")+ex.what()); 977 978 } 978 979 … … 992 993 // be adjacent, they are joined into the new one 993 994 void LFLineListOperations::addNewSearchResult(const std::list<pair<int, int> > &newlines, 994 std::list<std::pair<int, int> > &lines_list) 995 std::list<std::pair<int, int> > &lines_list) 995 996 throw(AipsError) 996 997 { … … 998 999 for (std::list<pair<int,int> >::const_iterator cli=newlines.begin(); 999 1000 cli!=newlines.end();++cli) { 1000 1001 1001 1002 // the first item, which has a non-void intersection or touches 1002 1003 // the new line 1003 1004 std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(), 1004 lines_list.end(), IntersectsWith(*cli)); 1005 // the last such item 1005 lines_list.end(), IntersectsWith(*cli)); 1006 // the last such item 1006 1007 std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg, 1007 1008 lines_list.end(), not1(IntersectsWith(*cli))); … … 1014 1015 lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end); 1015 1016 1016 // build a union of all intersecting lines 1017 // build a union of all intersecting lines 1017 1018 pair<int,int> union_line=for_each(lines_buffer.begin(), 1018 1019 lines_buffer.end(),BuildUnion(*cli)).result(); 1019 1020 1020 1021 // search for a right place for the new line (union_line) and add 1021 1022 std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(), … … 1026 1027 catch (const AipsError &ae) { 1027 1028 throw; 1028 } 1029 } 1029 1030 catch (const exception &ex) { 1030 1031 throw AipsError(String("LFLineListOperations::addNewSearchResult - STL error: ")+ex.what()); … … 1050 1051 if (signs[n]==signs[li->first] && signs[li->first]) 1051 1052 li->first=n; 1052 else break; 1053 else break; 1053 1054 } 1054 1055 // update the right hand side … … 1057 1058 if (signs[n]==signs[li->second-1] && signs[li->second-1]) 1058 1059 li->second=n; 1059 else break; 1060 else break; 1060 1061 } 1061 1062 } … … 1068 1069 catch (const AipsError &ae) { 1069 1070 throw; 1070 } 1071 } 1071 1072 catch (const exception &ex) { 1072 1073 throw AipsError(String("LFLineListOperations::extendLines - STL error: ")+ex.what()); -
trunk/src/SDLineFinder.h
r370 r881 1 1 //#--------------------------------------------------------------------------- 2 //# S DLineFinder.h: A class for automated spectral line search2 //# STLineFinder.h: A class for automated spectral line search 3 3 //#--------------------------------------------------------------------------- 4 4 //# Copyright (C) 2004 … … 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id: 29 //# $Id:$ 30 30 //#--------------------------------------------------------------------------- 31 #ifndef S DLINEFINDER_H32 #define S DLINEFINDER_H31 #ifndef STLINEFINDER_H 32 #define STLINEFINDER_H 33 33 34 34 // STL … … 49 49 50 50 // ASAP 51 #include "S DMemTableWrapper.h"52 #include "S DMemTable.h"51 #include "ScantableWrapper.h" 52 #include "Scantable.h" 53 53 54 54 namespace asap { … … 60 60 // 61 61 62 struct 62 struct LFLineListOperations { 63 63 // concatenate two lists preserving the order. If two lines appear to 64 // be adjacent or have a non-void intersection, they are joined into 64 // be adjacent or have a non-void intersection, they are joined into 65 65 // the new line 66 66 static void addNewSearchResult(const std::list<std::pair<int, int> > … … 79 79 throw(casa::AipsError); 80 80 protected: 81 81 82 82 // An auxiliary object function to test whether two lines have a non-void 83 83 // intersection … … 105 105 const std::pair<int,int>& result() const throw(); 106 106 }; 107 107 108 108 // An auxiliary object function to test whether a specified line 109 109 // is at lower spectral channels (to preserve the order in the line list) … … 117 117 // in the ordered list (so, it is at greater channel numbers) 118 118 bool operator()(const std::pair<int,int> &line2) const throw(); 119 }; 120 121 119 }; 120 121 122 122 }; 123 123 … … 127 127 /////////////////////////////////////////////////////////////////////////////// 128 128 // 129 // S DLineFinder - a class for automated spectral line search130 // 131 // 132 133 struct S DLineFinder : protected LFLineListOperations {134 S DLineFinder() throw();135 virtual ~S DLineFinder() throw(casa::AipsError);129 // STLineFinder - a class for automated spectral line search 130 // 131 // 132 133 struct STLineFinder : protected LFLineListOperations { 134 STLineFinder() throw(); 135 virtual ~STLineFinder() throw(casa::AipsError); 136 136 137 137 // set the parameters controlling algorithm … … 145 145 // in_avg_limit perform the averaging of no more than in_avg_limit 146 146 // adjacent channels to search for broad lines 147 // Default is 8, but for a bad baseline shape this 147 // Default is 8, but for a bad baseline shape this 148 148 // parameter should be decreased (may be even down to a 149 149 // minimum of 1 to disable this option) to avoid 150 150 // confusing of baseline undulations with a real line. 151 // Setting a very large value doesn't usually provide 152 // valid detections. 151 // Setting a very large value doesn't usually provide 152 // valid detections. 153 153 // in_box_size the box size for running mean calculation. Default is 154 154 // 1./5. of the whole spectrum size … … 164 164 // channels to drop from both sides of the spectrum 165 165 // in_edge is introduced for convinience, although all functionality 166 // can be achieved using a spectrum mask only 167 void setScan(const S DMemTableWrapper &in_scan,166 // can be achieved using a spectrum mask only 167 void setScan(const ScantableWrapper &in_scan, 168 168 const std::vector<bool> &in_mask, 169 169 const boost::python::tuple &in_edge) throw(casa::AipsError); … … 171 171 // search for spectral lines for a row specified by whichRow and 172 172 // Beam/IF/Pol specified by current cursor set for the scantable 173 // Number of lines found is returned 173 // Number of lines found is returned 174 174 int findLines(const casa::uInt &whichRow = 0) throw(casa::AipsError); 175 175 … … 182 182 183 183 // get range for all lines found. The same units as used in the scan 184 // will be returned (e.g. velocity instead of channels). 184 // will be returned (e.g. velocity instead of channels). 185 185 std::vector<double> getLineRanges() const throw(casa::AipsError); 186 186 // The same as getLineRanges, but channels are always used to specify … … 203 203 void subtractBaseline(const casa::Vector<casa::Bool> &temp_mask, 204 204 const casa::Int &order) throw(casa::AipsError); 205 205 206 206 // an auxiliary function to remove all lines from the list, except the 207 207 // strongest one (by absolute value). If the lines removed are real, 208 // they will be find again at the next iteration. This approach 209 // increases the number of iterations required, but is able to remove 208 // they will be find again at the next iteration. This approach 209 // increases the number of iterations required, but is able to remove 210 210 // the sidelobes likely to occur near strong lines. 211 211 // Later a better criterion may be implemented, e.g. 212 212 // taking into consideration the brightness of different lines. Now 213 // use the simplest solution 213 // use the simplest solution 214 214 // temp_mask - mask to work with (may be different from original mask as 215 215 // the lines previously found may be masked) … … 222 222 throw (casa::AipsError); 223 223 private: 224 casa::CountedConstPtr<S DMemTable> scan; // the scan to work with224 casa::CountedConstPtr<Scantable> scan; // the scan to work with 225 225 casa::Vector<casa::Bool> mask; // associated mask 226 226 std::pair<int,int> edge; // start and stop+1 channels 227 227 // to work with 228 casa::Float threshold; // detection threshold - the 228 casa::Float threshold; // detection threshold - the 229 229 // minimal signal to noise ratio 230 230 casa::Double box_size; // size of the box for running … … 251 251 252 252 } // namespace asap 253 #endif // #ifndef S DLINEFINDER_H253 #endif // #ifndef STLINEFINDER_H -
trunk/src/python_SDLineFinder.cc
r370 r881 1 1 //#--------------------------------------------------------------------------- 2 //# python_S DLineFinder.cc: python exposure of C++ SDLineFinder class2 //# python_STLineFinder.cc: python exposure of C++ STLineFinder class 3 3 //#--------------------------------------------------------------------------- 4 4 //# Copyright (C) 2004 … … 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id: 29 //# $Id:$ 30 30 //#--------------------------------------------------------------------------- 31 31 #include <boost/python.hpp> … … 37 37 namespace asap { 38 38 namespace python { 39 void python_S DLineFinder() {40 class_<S DLineFinder>("linefinder")39 void python_STLineFinder() { 40 class_<STLineFinder>("linefinder") 41 41 .def( init <> () ) 42 .def("setoptions",&S DLineFinder::setOptions)43 .def("setscan",&S DLineFinder::setScan)44 .def("findlines",&S DLineFinder::findLines)45 .def("getmask",&S DLineFinder::getMask)46 .def("getlineranges",&S DLineFinder::getLineRanges)47 .def("getlinerangesinchannels",&S DLineFinder::getLineRangesInChannels)42 .def("setoptions",&STLineFinder::setOptions) 43 .def("setscan",&STLineFinder::setScan) 44 .def("findlines",&STLineFinder::findLines) 45 .def("getmask",&STLineFinder::getMask) 46 .def("getlineranges",&STLineFinder::getLineRanges) 47 .def("getlinerangesinchannels",&STLineFinder::getLineRangesInChannels) 48 48 ; 49 49 }; -
trunk/src/python_asap.cpp
r875 r881 61 61 asap::python::python_STMath(); 62 62 asap::python::python_SDFitter(); 63 asap::python::python_STLineFinder(); 63 64 /* 64 65 asap::python::python_SDWriter(); 65 66 asap::python::python_SDFitTable(); 66 asap::python::python_SDLineFinder();67 67 */ 68 68 asap::python::python_SDLog(); -
trunk/src/python_asap.h
r875 r881 42 42 void python_STMath(); 43 43 void python_SDFitter(); 44 void python_STLineFinder(); 44 45 /* 45 46 void python_SDWriter(); 46 void python_SDMath();47 47 void python_SDFitTable(); 48 void python_SDLineFinder();49 48 */ 50 49 void python_SDLog();
Note:
See TracChangeset
for help on using the changeset viewer.