source: trunk/src/STLineFinder.cpp@ 1528

Last change on this file since 1528 was 1497, checked in by Max Voronkov, 16 years ago

fixed #149. The problem was inside getMask method. In some rare circumstances largely when single channel detections are allowed and with spectral averaging enabled, the method failed to flag the lines which have already been detected and went into infinite loop by detecting them again. The condition of leaving the loop (no new lines found) was never satisfied as a result of this.

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