source: trunk/src/SDLineFinder.cc@ 542

Last change on this file since 542 was 516, checked in by vor010, 20 years ago

LineFinder/automatic baseline fitter: a bug related to multiple-row scantable handling has been corrected.
Help is changed to describe a new interface. Parameter checking + vector to tuple conversion for the edge parameter have been added

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 38.8 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDLineFinder.cc: A class for automated spectral line search
3//#--------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id:
30//#---------------------------------------------------------------------------
31
32
33// ASAP
34#include "SDLineFinder.h"
35#include "SDFitter.h"
36
37// STL
38#include <functional>
39#include <algorithm>
40#include <iostream>
41#include <fstream>
42
43using namespace asap;
44using namespace casa;
45using namespace std;
46using namespace boost::python;
47
48namespace asap {
49
50///////////////////////////////////////////////////////////////////////////////
51//
52// RunningBox - a running box calculator. This class implements
53// interations over the specified spectrum and calculates
54// running box filter statistics.
55//
56
57class RunningBox {
58 // The input data to work with. Use reference symantics to avoid
59 // an unnecessary copying
60 const casa::Vector<casa::Float> &spectrum; // a buffer for the spectrum
61 const casa::Vector<casa::Bool> &mask; // associated mask
62 const std::pair<int,int> &edge; // start and stop+1 channels
63 // to work with
64
65 // statistics for running box filtering
66 casa::Float sumf; // sum of fluxes
67 casa::Float sumf2; // sum of squares of fluxes
68 casa::Float sumch; // sum of channel numbers (for linear fit)
69 casa::Float sumch2; // sum of squares of channel numbers (for linear fit)
70 casa::Float sumfch; // sum of flux*(channel number) (for linear fit)
71
72 int box_chan_cntr; // actual number of channels in the box
73 int max_box_nchan; // maximum allowed number of channels in the box
74 // (calculated from boxsize and actual spectrum size)
75 // cache for derivative statistics
76 mutable casa::Bool need2recalculate; // if true, values of the statistics
77 // below are invalid
78 mutable casa::Float linmean; // a value of the linear fit to the
79 // points in the running box
80 mutable casa::Float linvariance; // the same for variance
81 int cur_channel; // the number of the current channel
82 int start_advance; // number of channel from which the box can
83 // be moved (the middle of the box, if there is no
84 // masking)
85public:
86 // set up the object with the references to actual data
87 // as well as the number of channels in the running box
88 RunningBox(const casa::Vector<casa::Float> &in_spectrum,
89 const casa::Vector<casa::Bool> &in_mask,
90 const std::pair<int,int> &in_edge,
91 int in_max_box_nchan) throw(AipsError);
92
93 // access to the statistics
94 const casa::Float& getLinMean() const throw(AipsError);
95 const casa::Float& getLinVariance() const throw(AipsError);
96 const casa::Float aboveMean() const throw(AipsError);
97 int getChannel() const throw();
98
99 // actual number of channels in the box (max_box_nchan, if no channels
100 // are masked)
101 int getNumberOfBoxPoints() const throw();
102
103 // next channel
104 void next() throw(AipsError);
105
106 // checking whether there are still elements
107 casa::Bool haveMore() const throw();
108
109 // go to start
110 void rewind() throw(AipsError);
111
112protected:
113 // supplementary function to control running mean calculations.
114 // It adds a specified channel to the running mean box and
115 // removes (ch-maxboxnchan+1)'th channel from there
116 // Channels, for which the mask is false or index is beyond the
117 // allowed range, are ignored
118 void advanceRunningBox(int ch) throw(casa::AipsError);
119
120 // calculate derivative statistics. This function is const, because
121 // it updates the cache only
122 void updateDerivativeStatistics() const throw(AipsError);
123};
124
125//
126///////////////////////////////////////////////////////////////////////////////
127
128///////////////////////////////////////////////////////////////////////////////
129//
130// LFAboveThreshold An algorithm for line detection using running box
131// statistics. Line is detected if it is above the
132// specified threshold at the specified number of
133// consequtive channels. Prefix LF stands for Line Finder
134//
135class LFAboveThreshold : protected LFLineListOperations {
136 // temporary line edge channels and flag, which is True if the line
137 // was detected in the previous channels.
138 std::pair<int,int> cur_line;
139 casa::Bool is_detected_before;
140 int min_nchan; // A minimum number of consequtive
141 // channels, which should satisfy
142 // the detection criterion, to be
143 // a detection
144 casa::Float threshold; // detection threshold - the
145 // minimal signal to noise ratio
146 std::list<pair<int,int> > &lines; // list where detections are saved
147 // (pair: start and stop+1 channel)
148 RunningBox *running_box; // running box filter
149public:
150
151 // set up the detection criterion
152 LFAboveThreshold(std::list<pair<int,int> > &in_lines,
153 int in_min_nchan = 3,
154 casa::Float in_threshold = 5) throw();
155 virtual ~LFAboveThreshold() throw();
156
157 // replace the detection criterion
158 void setCriterion(int in_min_nchan, casa::Float in_threshold) throw();
159
160 // find spectral lines and add them into list
161 // if statholder is not NULL, the accumulate function of it will be
162 // called for each channel to save statistics
163 // spectrum, mask and edge - reference to the data
164 // max_box_nchan - number of channels in the running box
165 void findLines(const casa::Vector<casa::Float> &spectrum,
166 const casa::Vector<casa::Bool> &mask,
167 const std::pair<int,int> &edge,
168 int max_box_nchan) throw(casa::AipsError);
169
170protected:
171
172 // process a channel: update curline and is_detected before and
173 // add a new line to the list, if necessary using processCurLine()
174 // detect=true indicates that the current channel satisfies the criterion
175 void processChannel(Bool detect, const casa::Vector<casa::Bool> &mask)
176 throw(casa::AipsError);
177
178 // process the interval of channels stored in curline
179 // if it satisfies the criterion, add this interval as a new line
180 void processCurLine(const casa::Vector<casa::Bool> &mask)
181 throw(casa::AipsError);
182};
183
184//
185///////////////////////////////////////////////////////////////////////////////
186
187} // namespace asap
188
189///////////////////////////////////////////////////////////////////////////////
190//
191// RunningBox - a running box calculator. This class implements
192// interations over the specified spectrum and calculates
193// running box filter statistics.
194//
195
196// set up the object with the references to actual data
197// and the number of channels in the running box
198RunningBox::RunningBox(const casa::Vector<casa::Float> &in_spectrum,
199 const casa::Vector<casa::Bool> &in_mask,
200 const std::pair<int,int> &in_edge,
201 int in_max_box_nchan) throw(AipsError) :
202 spectrum(in_spectrum), mask(in_mask), edge(in_edge),
203 max_box_nchan(in_max_box_nchan)
204{
205 rewind();
206}
207
208void RunningBox::rewind() throw(AipsError) {
209 // fill statistics for initial box
210 box_chan_cntr=0; // no channels are currently in the box
211 sumf=0.; // initialize statistics
212 sumf2=0.;
213 sumch=0.;
214 sumch2=0.;
215 sumfch=0.;
216 int initial_box_ch=edge.first;
217 for (;initial_box_ch<edge.second && box_chan_cntr<max_box_nchan;
218 ++initial_box_ch)
219 advanceRunningBox(initial_box_ch);
220
221 if (initial_box_ch==edge.second)
222 throw AipsError("RunningBox::rewind - too much channels are masked");
223
224 cur_channel=edge.first;
225 start_advance=initial_box_ch-max_box_nchan/2;
226}
227
228// access to the statistics
229const casa::Float& RunningBox::getLinMean() const throw(AipsError)
230{
231 DebugAssert(cur_channel<edge.second, AipsError);
232 if (need2recalculate) updateDerivativeStatistics();
233 return linmean;
234}
235
236const casa::Float& RunningBox::getLinVariance() const throw(AipsError)
237{
238 DebugAssert(cur_channel<edge.second, AipsError);
239 if (need2recalculate) updateDerivativeStatistics();
240 return linvariance;
241}
242
243const casa::Float RunningBox::aboveMean() const throw(AipsError)
244{
245 DebugAssert(cur_channel<edge.second, AipsError);
246 if (need2recalculate) updateDerivativeStatistics();
247 return spectrum[cur_channel]-linmean;
248}
249
250int RunningBox::getChannel() const throw()
251{
252 return cur_channel;
253}
254
255// actual number of channels in the box (max_box_nchan, if no channels
256// are masked)
257int RunningBox::getNumberOfBoxPoints() const throw()
258{
259 return box_chan_cntr;
260}
261
262// supplementary function to control running mean calculations.
263// It adds a specified channel to the running mean box and
264// removes (ch-max_box_nchan+1)'th channel from there
265// Channels, for which the mask is false or index is beyond the
266// allowed range, are ignored
267void RunningBox::advanceRunningBox(int ch) throw(AipsError)
268{
269 if (ch>=edge.first && ch<edge.second)
270 if (mask[ch]) { // ch is a valid channel
271 ++box_chan_cntr;
272 sumf+=spectrum[ch];
273 sumf2+=square(spectrum[ch]);
274 sumch+=Float(ch);
275 sumch2+=square(Float(ch));
276 sumfch+=spectrum[ch]*Float(ch);
277 need2recalculate=True;
278 }
279 int ch2remove=ch-max_box_nchan;
280 if (ch2remove>=edge.first && ch2remove<edge.second)
281 if (mask[ch2remove]) { // ch2remove is a valid channel
282 --box_chan_cntr;
283 sumf-=spectrum[ch2remove];
284 sumf2-=square(spectrum[ch2remove]);
285 sumch-=Float(ch2remove);
286 sumch2-=square(Float(ch2remove));
287 sumfch-=spectrum[ch2remove]*Float(ch2remove);
288 need2recalculate=True;
289 }
290}
291
292// next channel
293void RunningBox::next() throw(AipsError)
294{
295 AlwaysAssert(cur_channel<edge.second,AipsError);
296 ++cur_channel;
297 if (cur_channel+max_box_nchan/2<edge.second && cur_channel>=start_advance)
298 advanceRunningBox(cur_channel+max_box_nchan/2); // update statistics
299}
300
301// checking whether there are still elements
302casa::Bool RunningBox::haveMore() const throw()
303{
304 return cur_channel<edge.second;
305}
306
307// calculate derivative statistics. This function is const, because
308// it updates the cache only
309void RunningBox::updateDerivativeStatistics() const throw(AipsError)
310{
311 AlwaysAssert(box_chan_cntr, AipsError);
312
313 Float mean=sumf/Float(box_chan_cntr);
314
315 // linear LSF formulae
316 Float meanch=sumch/Float(box_chan_cntr);
317 Float meanch2=sumch2/Float(box_chan_cntr);
318 if (meanch==meanch2 || box_chan_cntr<3) {
319 // vertical line in the spectrum, can't calculate linmean and linvariance
320 linmean=0.;
321 linvariance=0.;
322 } else {
323 Float coeff=(sumfch/Float(box_chan_cntr)-meanch*mean)/
324 (meanch2-square(meanch));
325 linmean=coeff*(Float(cur_channel)-meanch)+mean;
326 linvariance=sqrt(sumf2/Float(box_chan_cntr)-square(mean)-
327 square(coeff)*(meanch2-square(meanch)));
328 }
329 need2recalculate=False;
330}
331
332
333//
334///////////////////////////////////////////////////////////////////////////////
335
336///////////////////////////////////////////////////////////////////////////////
337//
338// LFAboveThreshold - a running mean algorithm for line detection
339//
340//
341
342
343// set up the detection criterion
344LFAboveThreshold::LFAboveThreshold(std::list<pair<int,int> > &in_lines,
345 int in_min_nchan,
346 casa::Float in_threshold) throw() :
347 min_nchan(in_min_nchan), threshold(in_threshold),
348 lines(in_lines), running_box(NULL) {}
349
350LFAboveThreshold::~LFAboveThreshold() throw()
351{
352 if (running_box!=NULL) delete running_box;
353}
354
355// replace the detection criterion
356void LFAboveThreshold::setCriterion(int in_min_nchan, casa::Float in_threshold)
357 throw()
358{
359 min_nchan=in_min_nchan;
360 threshold=in_threshold;
361}
362
363
364// process a channel: update cur_line and is_detected before and
365// add a new line to the list, if necessary
366void LFAboveThreshold::processChannel(Bool detect,
367 const casa::Vector<casa::Bool> &mask) throw(casa::AipsError)
368{
369 try {
370 if (detect) {
371 if (is_detected_before)
372 cur_line.second=running_box->getChannel()+1;
373 else {
374 is_detected_before=True;
375 cur_line.first=running_box->getChannel();
376 cur_line.second=running_box->getChannel()+1;
377 }
378 } else processCurLine(mask);
379 }
380 catch (const AipsError &ae) {
381 throw;
382 }
383 catch (const exception &ex) {
384 throw AipsError(String("LFAboveThreshold::processChannel - STL error: ")+ex.what());
385 }
386}
387
388// process the interval of channels stored in cur_line
389// if it satisfies the criterion, add this interval as a new line
390void LFAboveThreshold::processCurLine(const casa::Vector<casa::Bool> &mask)
391 throw(casa::AipsError)
392{
393 try {
394 if (is_detected_before) {
395 if (cur_line.second-cur_line.first>min_nchan) {
396 // it was a detection. We need to change the list
397 Bool add_new_line=False;
398 if (lines.size()) {
399 for (int i=lines.back().second;i<cur_line.first;++i)
400 if (mask[i]) { // one valid channel in between
401 // means that we deal with a separate line
402 add_new_line=True;
403 break;
404 }
405 } else add_new_line=True;
406 if (add_new_line)
407 lines.push_back(cur_line);
408 else lines.back().second=cur_line.second;
409 }
410 is_detected_before=False;
411 }
412 }
413 catch (const AipsError &ae) {
414 throw;
415 }
416 catch (const exception &ex) {
417 throw AipsError(String("LFAboveThreshold::processCurLine - STL error: ")+ex.what());
418 }
419}
420
421// find spectral lines and add them into list
422void LFAboveThreshold::findLines(const casa::Vector<casa::Float> &spectrum,
423 const casa::Vector<casa::Bool> &mask,
424 const std::pair<int,int> &edge,
425 int max_box_nchan)
426 throw(casa::AipsError)
427{
428 const int minboxnchan=4;
429 try {
430
431 if (running_box!=NULL) delete running_box;
432 running_box=new RunningBox(spectrum,mask,edge,max_box_nchan);
433
434
435 // determine the off-line variance first
436 // an assumption made: lines occupy a small part of the spectrum
437
438 std::vector<float> variances(edge.second-edge.first);
439 DebugAssert(variances.size(),AipsError);
440
441 for (;running_box->haveMore();running_box->next())
442 variances[running_box->getChannel()-edge.first]=
443 running_box->getLinVariance();
444
445 // in the future we probably should do a proper Chi^2 estimation
446 // now a simple 80% of smaller values will be used.
447 // it may degrade the performance of the algorithm for weak lines
448 // due to a bias of the Chi^2 distribution.
449 stable_sort(variances.begin(),variances.end());
450
451 Float offline_variance=0;
452 uInt offline_cnt=uInt(0.8*variances.size());
453 if (!offline_cnt) offline_cnt=variances.size(); // no much else left,
454 // although it is very inaccurate
455 for (uInt n=0;n<offline_cnt;++n)
456 offline_variance+=variances[n];
457 offline_variance/=Float(offline_cnt);
458
459 // actual search algorithm
460 is_detected_before=False;
461 Vector<Int> signs(spectrum.nelements(),0);
462
463 //ofstream os("dbg.dat");
464 for (running_box->rewind();running_box->haveMore();
465 running_box->next()) {
466 const int ch=running_box->getChannel();
467 if (running_box->getNumberOfBoxPoints()>=minboxnchan)
468 processChannel(mask[ch] && (fabs(running_box->aboveMean()) >=
469 threshold*offline_variance), mask);
470 else processCurLine(mask); // just finish what was accumulated before
471 const Float buf=running_box->aboveMean();
472 if (buf>0) signs[ch]=1;
473 else if (buf<0) signs[ch]=-1;
474 else if (buf==0) signs[ch]=0;
475 // os<<ch<<" "<<spectrum[ch]<<" "<<running_box->getLinMean()<<" "<<
476 // threshold*offline_variance<<endl;
477 }
478 if (lines.size())
479 searchForWings(lines,signs,mask,edge);
480 }
481 catch (const AipsError &ae) {
482 throw;
483 }
484 catch (const exception &ex) {
485 throw AipsError(String("LFAboveThreshold::findLines - STL error: ")+ex.what());
486 }
487}
488
489//
490///////////////////////////////////////////////////////////////////////////////
491
492///////////////////////////////////////////////////////////////////////////////
493//
494// LFLineListOperations::IntersectsWith - An auxiliary object function
495// to test whether two lines have a non-void intersection
496//
497
498
499// line1 - range of the first line: start channel and stop+1
500LFLineListOperations::IntersectsWith::IntersectsWith(const std::pair<int,int> &in_line1) :
501 line1(in_line1) {}
502
503
504// return true if line2 intersects with line1 with at least one
505// common channel, and false otherwise
506// line2 - range of the second line: start channel and stop+1
507bool LFLineListOperations::IntersectsWith::operator()(const std::pair<int,int> &line2)
508 const throw()
509{
510 if (line2.second<line1.first) return false; // line2 is at lower channels
511 if (line2.first>line1.second) return false; // line2 is at upper channels
512 return true; // line2 has an intersection or is adjacent to line1
513}
514
515//
516///////////////////////////////////////////////////////////////////////////////
517
518///////////////////////////////////////////////////////////////////////////////
519//
520// LFLineListOperations::BuildUnion - An auxiliary object function to build a union
521// of several lines to account for a possibility of merging the nearby lines
522//
523
524// set an initial line (can be a first line in the sequence)
525LFLineListOperations::BuildUnion::BuildUnion(const std::pair<int,int> &line1) :
526 temp_line(line1) {}
527
528// update temp_line with a union of temp_line and new_line
529// provided there is no gap between the lines
530void LFLineListOperations::BuildUnion::operator()(const std::pair<int,int> &new_line)
531 throw()
532{
533 if (new_line.first<temp_line.first) temp_line.first=new_line.first;
534 if (new_line.second>temp_line.second) temp_line.second=new_line.second;
535}
536
537// return the result (temp_line)
538const std::pair<int,int>& LFLineListOperations::BuildUnion::result() const throw()
539{
540 return temp_line;
541}
542
543//
544///////////////////////////////////////////////////////////////////////////////
545
546///////////////////////////////////////////////////////////////////////////////
547//
548// LFLineListOperations::LaterThan - An auxiliary object function to test whether a
549// specified line is at lower spectral channels (to preserve the order in
550// the line list)
551//
552
553// setup the line to compare with
554LFLineListOperations::LaterThan::LaterThan(const std::pair<int,int> &in_line1) :
555 line1(in_line1) {}
556
557// return true if line2 should be placed later than line1
558// in the ordered list (so, it is at greater channel numbers)
559bool LFLineListOperations::LaterThan::operator()(const std::pair<int,int> &line2)
560 const throw()
561{
562 if (line2.second<line1.first) return false; // line2 is at lower channels
563 if (line2.first>line1.second) return true; // line2 is at upper channels
564
565 // line2 intersects with line1. We should have no such situation in
566 // practice
567 return line2.first>line1.first;
568}
569
570//
571///////////////////////////////////////////////////////////////////////////////
572
573
574///////////////////////////////////////////////////////////////////////////////
575//
576// SDLineFinder - a class for automated spectral line search
577//
578//
579
580SDLineFinder::SDLineFinder() throw() : edge(0,0)
581{
582 setOptions();
583}
584
585// set the parameters controlling algorithm
586// in_threshold a single channel threshold default is sqrt(3), which
587// means together with 3 minimum channels at least 3 sigma
588// detection criterion
589// For bad baseline shape, in_threshold may need to be
590// increased
591// in_min_nchan minimum number of channels above the threshold to report
592// a detection, default is 3
593// in_avg_limit perform the averaging of no more than in_avg_limit
594// adjacent channels to search for broad lines
595// Default is 8, but for a bad baseline shape this
596// parameter should be decreased (may be even down to a
597// minimum of 1 to disable this option) to avoid
598// confusing of baseline undulations with a real line.
599// Setting a very large value doesn't usually provide
600// valid detections.
601// in_box_size the box size for running mean calculation. Default is
602// 1./5. of the whole spectrum size
603void SDLineFinder::setOptions(const casa::Float &in_threshold,
604 const casa::Int &in_min_nchan,
605 const casa::Int &in_avg_limit,
606 const casa::Float &in_box_size) throw()
607{
608 threshold=in_threshold;
609 min_nchan=in_min_nchan;
610 avg_limit=in_avg_limit;
611 box_size=in_box_size;
612}
613
614SDLineFinder::~SDLineFinder() throw(AipsError) {}
615
616// set scan to work with (in_scan parameter), associated mask (in_mask
617// parameter) and the edge channel rejection (in_edge parameter)
618// if in_edge has zero length, all channels chosen by mask will be used
619// if in_edge has one element only, it represents the number of
620// channels to drop from both sides of the spectrum
621// in_edge is introduced for convinience, although all functionality
622// can be achieved using a spectrum mask only
623void SDLineFinder::setScan(const SDMemTableWrapper &in_scan,
624 const std::vector<bool> &in_mask,
625 const boost::python::tuple &in_edge) throw(AipsError)
626{
627 try {
628 scan=in_scan.getCP();
629 AlwaysAssert(!scan.null(),AipsError);
630
631 mask=in_mask;
632 if (mask.nelements()!=scan->nChan())
633 throw AipsError("SDLineFinder::setScan - in_scan and in_mask have different"
634 "number of spectral channels.");
635
636 // number of elements in the in_edge tuple
637 int n=extract<int>(in_edge.attr("__len__")());
638 if (n>2 || n<0)
639 throw AipsError("SDLineFinder::setScan - the length of the in_edge parameter"
640 "should not exceed 2");
641 if (!n) {
642 // all spectrum, no rejection
643 edge.first=0;
644 edge.second=scan->nChan();
645 } else {
646 edge.first=extract<int>(in_edge.attr("__getitem__")(0));
647 if (edge.first<0)
648 throw AipsError("SDLineFinder::setScan - the in_edge parameter has a negative"
649 "number of channels to drop");
650 if (edge.first>=scan->nChan())
651 throw AipsError("SDLineFinder::setScan - all channels are rejected by the in_edge parameter");
652 if (n==2) {
653 edge.second=extract<int>(in_edge.attr("__getitem__")(1));
654 if (edge.second<0)
655 throw AipsError("SDLineFinder::setScan - the in_edge parameter has a negative"
656 "number of channels to drop");
657 edge.second=scan->nChan()-edge.second;
658 } else edge.second=scan->nChan()-edge.first;
659 if (edge.second<0 || (edge.first>=edge.second))
660 throw AipsError("SDLineFinder::setScan - all channels are rejected by the in_edge parameter");
661 }
662 }
663 catch(const AipsError &ae) {
664 // setScan is unsuccessfull, reset scan/mask/edge
665 scan=CountedConstPtr<SDMemTable>(); // null pointer
666 mask.resize(0);
667 edge=pair<int,int>(0,0);
668 throw;
669 }
670}
671
672// search for spectral lines. Number of lines found is returned
673int SDLineFinder::findLines(const casa::uInt &whichRow) throw(casa::AipsError)
674{
675 const int minboxnchan=4;
676 if (scan.null())
677 throw AipsError("SDLineFinder::findLines - a scan should be set first,"
678 " use set_scan");
679 DebugAssert(mask.nelements()==scan->nChan(), AipsError);
680 int max_box_nchan=int(scan->nChan()*box_size); // number of channels in running
681 // box
682 if (max_box_nchan<2)
683 throw AipsError("SDLineFinder::findLines - box_size is too small");
684
685 scan->getSpectrum(spectrum, whichRow);
686
687 lines.resize(0); // search from the scratch
688 last_row_used=whichRow;
689 Vector<Bool> temp_mask(mask);
690
691 Bool first_pass=True;
692 Int avg_factor=1; // this number of adjacent channels is averaged together
693 // the total number of the channels is not altered
694 // instead, min_nchan is also scaled
695 // it helps to search for broad lines
696 while (true) {
697 // a buffer for new lines found at this iteration
698 std::list<pair<int,int> > new_lines;
699
700 try {
701 // line find algorithm
702 LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold);
703 lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan);
704 first_pass=False;
705 if (!new_lines.size())
706 throw AipsError("spurious"); // nothing new - use the same
707 // code as for a real exception
708 }
709 catch(const AipsError &ae) {
710 if (first_pass) throw;
711 // nothing new - proceed to the next step of averaging, if any
712 // (to search for broad lines)
713 if (avg_factor>avg_limit) break; // averaging up to avg_limit
714 // adjacent channels,
715 // stop after that
716 avg_factor*=2; // twice as more averaging
717 subtractBaseline(temp_mask,9);
718 averageAdjacentChannels(temp_mask,avg_factor);
719 continue;
720 }
721 keepStrongestOnly(temp_mask,new_lines,max_box_nchan);
722 // update the list (lines) merging intervals, if necessary
723 addNewSearchResult(new_lines,lines);
724 // get a new mask
725 temp_mask=getMask();
726 }
727 return int(lines.size());
728}
729
730// auxiliary function to fit and subtract a polynomial from the current
731// spectrum. It uses the SDFitter class. This action is required before
732// reducing the spectral resolution if the baseline shape is bad
733void SDLineFinder::subtractBaseline(const casa::Vector<casa::Bool> &temp_mask,
734 const casa::Int &order) throw(casa::AipsError)
735{
736 AlwaysAssert(spectrum.nelements(),AipsError);
737 // use the fact that temp_mask excludes channels rejected at the edge
738 SDFitter sdf;
739 std::vector<float> absc(spectrum.nelements());
740 for (Int i=0;i<absc.size();++i)
741 absc[i]=float(i)/float(spectrum.nelements());
742 std::vector<float> spec;
743 spectrum.tovector(spec);
744 std::vector<bool> std_mask;
745 temp_mask.tovector(std_mask);
746 sdf.setData(absc,spec,std_mask);
747 sdf.setExpression("poly",order);
748 if (!sdf.fit()) return; // fit failed, use old spectrum
749 spectrum=casa::Vector<casa::Float>(sdf.getResidual());
750}
751
752// auxiliary function to average adjacent channels and update the mask
753// if at least one channel involved in summation is masked, all
754// output channels will be masked. This function works with the
755// spectrum and edge fields of this class, but updates the mask
756// array specified, rather than the field of this class
757// boxsize - a number of adjacent channels to average
758void SDLineFinder::averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update,
759 const casa::Int &boxsize)
760 throw(casa::AipsError)
761{
762 DebugAssert(mask2update.nelements()==spectrum.nelements(), AipsError);
763 DebugAssert(boxsize!=0,AipsError);
764
765 for (int n=edge.first;n<edge.second;n+=boxsize) {
766 DebugAssert(n<spectrum.nelements(),AipsError);
767 int nboxch=0; // number of channels currently in the box
768 Float mean=0; // buffer for mean calculations
769 for (int k=n;k<n+boxsize && k<edge.second;++k)
770 if (mask2update[k]) { // k is a valid channel
771 mean+=spectrum[k];
772 ++nboxch;
773 }
774 if (nboxch<boxsize) // mask these channels
775 for (int k=n;k<n+boxsize && k<edge.second;++k)
776 mask2update[k]=False;
777 else {
778 mean/=Float(boxsize);
779 for (int k=n;k<n+boxsize && k<edge.second;++k)
780 spectrum[k]=mean;
781 }
782 }
783}
784
785
786// get the mask to mask out all lines that have been found (default)
787// if invert=true, only channels belong to lines will be unmasked
788// Note: all channels originally masked by the input mask (in_mask
789// in setScan) or dropped out by the edge parameter (in_edge
790// in setScan) are still excluded regardless on the invert option
791std::vector<bool> SDLineFinder::getMask(bool invert)
792 const throw(casa::AipsError)
793{
794 try {
795 if (scan.null())
796 throw AipsError("SDLineFinder::getMask - a scan should be set first,"
797 " use set_scan followed by find_lines");
798 DebugAssert(mask.nelements()==scan->nChan(), AipsError);
799 /*
800 if (!lines.size())
801 throw AipsError("SDLineFinder::getMask - one have to search for "
802 "lines first, use find_lines");
803 */
804 std::vector<bool> res_mask(mask.nelements());
805 // iterator through lines
806 std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
807 for (int ch=0;ch<res_mask.size();++ch)
808 if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
809 else if (!mask[ch]) res_mask[ch]=false;
810 else {
811 res_mask[ch]=!invert; // no line by default
812 if (cli==lines.end()) continue;
813 if (ch>=cli->first && ch<cli->second)
814 res_mask[ch]=invert; // this is a line
815 if (ch>=cli->second)
816 ++cli; // next line in the list
817 }
818
819 return res_mask;
820 }
821 catch (const AipsError &ae) {
822 throw;
823 }
824 catch (const exception &ex) {
825 throw AipsError(String("SDLineFinder::getMask - STL error: ")+ex.what());
826 }
827}
828
829// get range for all lines found. The same units as used in the scan
830// will be returned (e.g. velocity instead of channels).
831std::vector<double> SDLineFinder::getLineRanges()
832 const throw(casa::AipsError)
833{
834 // convert to required abscissa units
835 std::vector<double> vel=scan->getAbcissa(last_row_used);
836 std::vector<int> ranges=getLineRangesInChannels();
837 std::vector<double> res(ranges.size());
838
839 std::vector<int>::const_iterator cri=ranges.begin();
840 std::vector<double>::iterator outi=res.begin();
841 for (;cri!=ranges.end() && outi!=res.end();++cri,++outi)
842 if (uInt(*cri)>=vel.size())
843 throw AipsError("SDLineFinder::getLineRanges - getAbcissa provided less channels than reqired");
844 else *outi=vel[*cri];
845 return res;
846}
847
848// The same as getLineRanges, but channels are always used to specify
849// the range
850std::vector<int> SDLineFinder::getLineRangesInChannels()
851 const throw(casa::AipsError)
852{
853 try {
854 if (scan.null())
855 throw AipsError("SDLineFinder::getLineRangesInChannels - a scan should be set first,"
856 " use set_scan followed by find_lines");
857 DebugAssert(mask.nelements()==scan->nChan(), AipsError);
858
859 if (!lines.size())
860 throw AipsError("SDLineFinder::getLineRangesInChannels - one have to search for "
861 "lines first, use find_lines");
862
863 std::vector<int> res(2*lines.size());
864 // iterator through lines & result
865 std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
866 std::vector<int>::iterator ri=res.begin();
867 for (;cli!=lines.end() && ri!=res.end();++cli,++ri) {
868 *ri=cli->first;
869 if (++ri!=res.end())
870 *ri=cli->second-1;
871 }
872 return res;
873 }
874 catch (const AipsError &ae) {
875 throw;
876 }
877 catch (const exception &ex) {
878 throw AipsError(String("SDLineFinder::getLineRanges - STL error: ")+ex.what());
879 }
880}
881
882
883
884// an auxiliary function to remove all lines from the list, except the
885// strongest one (by absolute value). If the lines removed are real,
886// they will be find again at the next iteration. This approach
887// increases the number of iterations required, but is able to remove
888// the sidelobes likely to occur near strong lines.
889// Later a better criterion may be implemented, e.g.
890// taking into consideration the brightness of different lines. Now
891// use the simplest solution
892// temp_mask - mask to work with (may be different from original mask as
893// the lines previously found may be masked)
894// lines2update - a list of lines to work with
895// nothing will be done if it is empty
896// max_box_nchan - channels in the running box for baseline filtering
897void SDLineFinder::keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask,
898 std::list<std::pair<int, int> > &lines2update,
899 int max_box_nchan)
900 throw (casa::AipsError)
901{
902 try {
903 if (!lines2update.size()) return; // ignore an empty list
904
905 // current line
906 std::list<std::pair<int,int> >::iterator li=lines2update.begin();
907 // strongest line
908 std::list<std::pair<int,int> >::iterator strongli=lines2update.begin();
909 // the flux (absolute value) of the strongest line
910 Float peak_flux=-1; // negative value - a flag showing uninitialized
911 // value
912 // the algorithm below relies on the list being ordered
913 Float tmp_flux=-1; // a temporary peak
914 for (RunningBox running_box(spectrum,temp_mask,edge,max_box_nchan);
915 running_box.haveMore(); running_box.next()) {
916
917 if (li==lines2update.end()) break; // no more lines
918 const int ch=running_box.getChannel();
919 if (ch>=li->first && ch<li->second)
920 if (temp_mask[ch] && tmp_flux<fabs(running_box.aboveMean()))
921 tmp_flux=fabs(running_box.aboveMean());
922 if (ch==li->second-1) {
923 if (peak_flux<tmp_flux) { // if peak_flux=-1, this condition
924 peak_flux=tmp_flux; // will be satisfied
925 strongli=li;
926 }
927 ++li;
928 tmp_flux=-1;
929 }
930 }
931 std::list<std::pair<int,int> > res;
932 res.splice(res.end(),lines2update,strongli);
933 lines2update.clear();
934 lines2update.splice(lines2update.end(),res);
935 }
936 catch (const AipsError &ae) {
937 throw;
938 }
939 catch (const exception &ex) {
940 throw AipsError(String("SDLineFinder::keepStrongestOnly - STL error: ")+ex.what());
941 }
942
943}
944
945//
946///////////////////////////////////////////////////////////////////////////////
947
948
949///////////////////////////////////////////////////////////////////////////////
950//
951// LFLineListOperations - a class incapsulating operations with line lists
952// The LF prefix stands for Line Finder
953//
954
955// concatenate two lists preserving the order. If two lines appear to
956// be adjacent, they are joined into the new one
957void LFLineListOperations::addNewSearchResult(const std::list<pair<int, int> > &newlines,
958 std::list<std::pair<int, int> > &lines_list)
959 throw(AipsError)
960{
961 try {
962 for (std::list<pair<int,int> >::const_iterator cli=newlines.begin();
963 cli!=newlines.end();++cli) {
964
965 // the first item, which has a non-void intersection or touches
966 // the new line
967 std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(),
968 lines_list.end(), IntersectsWith(*cli));
969 // the last such item
970 std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg,
971 lines_list.end(), not1(IntersectsWith(*cli)));
972
973 // extract all lines which intersect or touch a new one into
974 // a temporary buffer. This may invalidate the iterators
975 // line_buffer may be empty, if no lines intersects with a new
976 // one.
977 std::list<pair<int,int> > lines_buffer;
978 lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end);
979
980 // build a union of all intersecting lines
981 pair<int,int> union_line=for_each(lines_buffer.begin(),
982 lines_buffer.end(),BuildUnion(*cli)).result();
983
984 // search for a right place for the new line (union_line) and add
985 std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(),
986 lines_list.end(), LaterThan(union_line));
987 lines_list.insert(pos2insert,union_line);
988 }
989 }
990 catch (const AipsError &ae) {
991 throw;
992 }
993 catch (const exception &ex) {
994 throw AipsError(String("LFLineListOperations::addNewSearchResult - STL error: ")+ex.what());
995 }
996}
997
998// extend all line ranges to the point where a value stored in the
999// specified vector changes (e.g. value-mean change its sign)
1000// This operation is necessary to include line wings, which are below
1001// the detection threshold. If lines becomes adjacent, they are
1002// merged together. Any masked channel stops the extension
1003void LFLineListOperations::searchForWings(std::list<std::pair<int, int> > &newlines,
1004 const casa::Vector<casa::Int> &signs,
1005 const casa::Vector<casa::Bool> &mask,
1006 const std::pair<int,int> &edge) throw(casa::AipsError)
1007{
1008 try {
1009 for (std::list<pair<int,int> >::iterator li=newlines.begin();
1010 li!=newlines.end();++li) {
1011 // update the left hand side
1012 for (int n=li->first-1;n>=edge.first;--n) {
1013 if (!mask[n]) break;
1014 if (signs[n]==signs[li->first] && signs[li->first])
1015 li->first=n;
1016 else break;
1017 }
1018 // update the right hand side
1019 for (int n=li->second;n<edge.second;++n) {
1020 if (!mask[n]) break;
1021 if (signs[n]==signs[li->second-1] && signs[li->second-1])
1022 li->second=n;
1023 else break;
1024 }
1025 }
1026 // need to search for possible mergers.
1027 std::list<std::pair<int, int> > result_buffer;
1028 addNewSearchResult(newlines,result_buffer);
1029 newlines.clear();
1030 newlines.splice(newlines.end(),result_buffer);
1031 }
1032 catch (const AipsError &ae) {
1033 throw;
1034 }
1035 catch (const exception &ex) {
1036 throw AipsError(String("LFLineListOperations::extendLines - STL error: ")+ex.what());
1037 }
1038}
1039
1040//
1041///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.