source: trunk/src/STLineFinder.cpp@ 965

Last change on this file since 965 was 924, checked in by mar637, 19 years ago

fix for incorrect use of Scantable::nchan. Fix for std::vector / casa::Vector conusion

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 40.9 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 924 2006-03-24 02:33:17Z mar637 $
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;
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
149 casa::Vector<Int> signs; // An array to store the signs of
150 // the value - current mean
151 // (used to search wings)
152 casa::Int last_sign; // a sign (+1, -1 or 0) of the
153 // last point of the detected line
154 //
155public:
156
157 // set up the detection criterion
158 LFAboveThreshold(std::list<pair<int,int> > &in_lines,
159 int in_min_nchan = 3,
160 casa::Float in_threshold = 5) throw();
161 virtual ~LFAboveThreshold() throw();
162
163 // replace the detection criterion
164 void setCriterion(int in_min_nchan, casa::Float in_threshold) throw();
165
166 // return the array with signs of the value-current mean
167 // An element is +1 if value>mean, -1 if less, 0 if equal.
168 // This array is updated each time the findLines method is called and
169 // is used to search the line wings
170 const casa::Vector<Int>& getSigns() const throw();
171
172 // find spectral lines and add them into list
173 // if statholder is not NULL, the accumulate function of it will be
174 // called for each channel to save statistics
175 // spectrum, mask and edge - reference to the data
176 // max_box_nchan - number of channels in the running box
177 void findLines(const casa::Vector<casa::Float> &spectrum,
178 const casa::Vector<casa::Bool> &mask,
179 const std::pair<int,int> &edge,
180 int max_box_nchan) throw(casa::AipsError);
181
182protected:
183
184 // process a channel: update curline and is_detected before and
185 // add a new line to the list, if necessary using processCurLine()
186 // detect=true indicates that the current channel satisfies the criterion
187 void processChannel(Bool detect, const casa::Vector<casa::Bool> &mask)
188 throw(casa::AipsError);
189
190 // process the interval of channels stored in curline
191 // if it satisfies the criterion, add this interval as a new line
192 void processCurLine(const casa::Vector<casa::Bool> &mask)
193 throw(casa::AipsError);
194
195 // get the sign of runningBox->aboveMean(). The RunningBox pointer
196 // should be defined
197 casa::Int getAboveMeanSign() const throw();
198};
199
200//
201///////////////////////////////////////////////////////////////////////////////
202
203} // namespace asap
204
205///////////////////////////////////////////////////////////////////////////////
206//
207// RunningBox - a running box calculator. This class implements
208// interations over the specified spectrum and calculates
209// running box filter statistics.
210//
211
212// set up the object with the references to actual data
213// and the number of channels in the running box
214RunningBox::RunningBox(const casa::Vector<casa::Float> &in_spectrum,
215 const casa::Vector<casa::Bool> &in_mask,
216 const std::pair<int,int> &in_edge,
217 int in_max_box_nchan) throw(AipsError) :
218 spectrum(in_spectrum), mask(in_mask), edge(in_edge),
219 max_box_nchan(in_max_box_nchan)
220{
221 rewind();
222}
223
224void RunningBox::rewind() throw(AipsError) {
225 // fill statistics for initial box
226 box_chan_cntr=0; // no channels are currently in the box
227 sumf=0.; // initialize statistics
228 sumf2=0.;
229 sumch=0.;
230 sumch2=0.;
231 sumfch=0.;
232 int initial_box_ch=edge.first;
233 for (;initial_box_ch<edge.second && box_chan_cntr<max_box_nchan;
234 ++initial_box_ch)
235 advanceRunningBox(initial_box_ch);
236
237 if (initial_box_ch==edge.second)
238 throw AipsError("RunningBox::rewind - too much channels are masked");
239
240 cur_channel=edge.first;
241 start_advance=initial_box_ch-max_box_nchan/2;
242}
243
244// access to the statistics
245const casa::Float& RunningBox::getLinMean() const throw(AipsError)
246{
247 DebugAssert(cur_channel<edge.second, AipsError);
248 if (need2recalculate) updateDerivativeStatistics();
249 return linmean;
250}
251
252const casa::Float& RunningBox::getLinVariance() const throw(AipsError)
253{
254 DebugAssert(cur_channel<edge.second, AipsError);
255 if (need2recalculate) updateDerivativeStatistics();
256 return linvariance;
257}
258
259const casa::Float RunningBox::aboveMean() const throw(AipsError)
260{
261 DebugAssert(cur_channel<edge.second, AipsError);
262 if (need2recalculate) updateDerivativeStatistics();
263 return spectrum[cur_channel]-linmean;
264}
265
266int RunningBox::getChannel() const throw()
267{
268 return cur_channel;
269}
270
271// actual number of channels in the box (max_box_nchan, if no channels
272// are masked)
273int RunningBox::getNumberOfBoxPoints() const throw()
274{
275 return box_chan_cntr;
276}
277
278// supplementary function to control running mean calculations.
279// It adds a specified channel to the running mean box and
280// removes (ch-max_box_nchan+1)'th channel from there
281// Channels, for which the mask is false or index is beyond the
282// allowed range, are ignored
283void RunningBox::advanceRunningBox(int ch) throw(AipsError)
284{
285 if (ch>=edge.first && ch<edge.second)
286 if (mask[ch]) { // ch is a valid channel
287 ++box_chan_cntr;
288 sumf+=spectrum[ch];
289 sumf2+=square(spectrum[ch]);
290 sumch+=Float(ch);
291 sumch2+=square(Float(ch));
292 sumfch+=spectrum[ch]*Float(ch);
293 need2recalculate=True;
294 }
295 int ch2remove=ch-max_box_nchan;
296 if (ch2remove>=edge.first && ch2remove<edge.second)
297 if (mask[ch2remove]) { // ch2remove is a valid channel
298 --box_chan_cntr;
299 sumf-=spectrum[ch2remove];
300 sumf2-=square(spectrum[ch2remove]);
301 sumch-=Float(ch2remove);
302 sumch2-=square(Float(ch2remove));
303 sumfch-=spectrum[ch2remove]*Float(ch2remove);
304 need2recalculate=True;
305 }
306}
307
308// next channel
309void RunningBox::next() throw(AipsError)
310{
311 AlwaysAssert(cur_channel<edge.second,AipsError);
312 ++cur_channel;
313 if (cur_channel+max_box_nchan/2<edge.second && cur_channel>=start_advance)
314 advanceRunningBox(cur_channel+max_box_nchan/2); // update statistics
315}
316
317// checking whether there are still elements
318casa::Bool RunningBox::haveMore() const throw()
319{
320 return cur_channel<edge.second;
321}
322
323// calculate derivative statistics. This function is const, because
324// it updates the cache only
325void RunningBox::updateDerivativeStatistics() const throw(AipsError)
326{
327 AlwaysAssert(box_chan_cntr, AipsError);
328
329 Float mean=sumf/Float(box_chan_cntr);
330
331 // linear LSF formulae
332 Float meanch=sumch/Float(box_chan_cntr);
333 Float meanch2=sumch2/Float(box_chan_cntr);
334 if (meanch==meanch2 || box_chan_cntr<3) {
335 // vertical line in the spectrum, can't calculate linmean and linvariance
336 linmean=0.;
337 linvariance=0.;
338 } else {
339 Float coeff=(sumfch/Float(box_chan_cntr)-meanch*mean)/
340 (meanch2-square(meanch));
341 linmean=coeff*(Float(cur_channel)-meanch)+mean;
342 linvariance=sqrt(sumf2/Float(box_chan_cntr)-square(mean)-
343 square(coeff)*(meanch2-square(meanch)));
344 }
345 need2recalculate=False;
346}
347
348
349//
350///////////////////////////////////////////////////////////////////////////////
351
352///////////////////////////////////////////////////////////////////////////////
353//
354// LFAboveThreshold - a running mean algorithm for line detection
355//
356//
357
358
359// set up the detection criterion
360LFAboveThreshold::LFAboveThreshold(std::list<pair<int,int> > &in_lines,
361 int in_min_nchan,
362 casa::Float in_threshold) throw() :
363 min_nchan(in_min_nchan), threshold(in_threshold),
364 lines(in_lines), running_box(NULL) {}
365
366LFAboveThreshold::~LFAboveThreshold() throw()
367{
368 if (running_box!=NULL) delete running_box;
369}
370
371// replace the detection criterion
372void LFAboveThreshold::setCriterion(int in_min_nchan, casa::Float in_threshold)
373 throw()
374{
375 min_nchan=in_min_nchan;
376 threshold=in_threshold;
377}
378
379// get the sign of runningBox->aboveMean(). The RunningBox pointer
380// should be defined
381casa::Int LFAboveThreshold::getAboveMeanSign() const throw()
382{
383 const Float buf=running_box->aboveMean();
384 if (buf>0) return 1;
385 if (buf<0) return -1;
386 return 0;
387}
388
389
390// process a channel: update cur_line and is_detected before and
391// add a new line to the list, if necessary
392void LFAboveThreshold::processChannel(Bool detect,
393 const casa::Vector<casa::Bool> &mask) throw(casa::AipsError)
394{
395 try {
396 if (is_detected_before) {
397 // we have to check that the current detection has the
398 // same sign of running_box->aboveMean
399 // otherwise it could be a spurious detection
400 if (last_sign && last_sign!=getAboveMeanSign())
401 detect=False;
402 }
403 if (detect) {
404 last_sign=getAboveMeanSign();
405 if (is_detected_before)
406 cur_line.second=running_box->getChannel()+1;
407 else {
408 is_detected_before=True;
409 cur_line.first=running_box->getChannel();
410 cur_line.second=running_box->getChannel()+1;
411 }
412 } else processCurLine(mask);
413 }
414 catch (const AipsError &ae) {
415 throw;
416 }
417 catch (const exception &ex) {
418 throw AipsError(String("LFAboveThreshold::processChannel - STL error: ")+ex.what());
419 }
420}
421
422// process the interval of channels stored in cur_line
423// if it satisfies the criterion, add this interval as a new line
424void LFAboveThreshold::processCurLine(const casa::Vector<casa::Bool> &mask)
425 throw(casa::AipsError)
426{
427 try {
428 if (is_detected_before) {
429 if (cur_line.second-cur_line.first>min_nchan) {
430 // it was a detection. We need to change the list
431 Bool add_new_line=False;
432 if (lines.size()) {
433 for (int i=lines.back().second;i<cur_line.first;++i)
434 if (mask[i]) { // one valid channel in between
435 // means that we deal with a separate line
436 add_new_line=True;
437 break;
438 }
439 } else add_new_line=True;
440 if (add_new_line)
441 lines.push_back(cur_line);
442 else lines.back().second=cur_line.second;
443 }
444 is_detected_before=False;
445 }
446 }
447 catch (const AipsError &ae) {
448 throw;
449 }
450 catch (const exception &ex) {
451 throw AipsError(String("LFAboveThreshold::processCurLine - STL error: ")+ex.what());
452 }
453}
454
455// return the array with signs of the value-current mean
456// An element is +1 if value>mean, -1 if less, 0 if equal.
457// This array is updated each time the findLines method is called and
458// is used to search the line wings
459const casa::Vector<Int>& LFAboveThreshold::getSigns() const throw()
460{
461 return signs;
462}
463
464// find spectral lines and add them into list
465void LFAboveThreshold::findLines(const casa::Vector<casa::Float> &spectrum,
466 const casa::Vector<casa::Bool> &mask,
467 const std::pair<int,int> &edge,
468 int max_box_nchan)
469 throw(casa::AipsError)
470{
471 const int minboxnchan=4;
472 try {
473
474 if (running_box!=NULL) delete running_box;
475 running_box=new RunningBox(spectrum,mask,edge,max_box_nchan);
476
477
478 // determine the off-line variance first
479 // an assumption made: lines occupy a small part of the spectrum
480
481 std::vector<float> variances(edge.second-edge.first);
482 DebugAssert(variances.size(),AipsError);
483
484 for (;running_box->haveMore();running_box->next())
485 variances[running_box->getChannel()-edge.first]=
486 running_box->getLinVariance();
487
488 // in the future we probably should do a proper Chi^2 estimation
489 // now a simple 80% of smaller values will be used.
490 // it may degrade the performance of the algorithm for weak lines
491 // due to a bias of the Chi^2 distribution.
492 stable_sort(variances.begin(),variances.end());
493
494 Float offline_variance=0;
495 uInt offline_cnt=uInt(0.8*variances.size());
496 if (!offline_cnt) offline_cnt=variances.size(); // no much else left,
497 // although it is very inaccurate
498 for (uInt n=0;n<offline_cnt;++n)
499 offline_variance+=variances[n];
500 offline_variance/=Float(offline_cnt);
501
502 // actual search algorithm
503 is_detected_before=False;
504
505 // initiate the signs array
506 signs.resize(spectrum.nelements());
507 signs=Vector<Int>(spectrum.nelements(),0);
508
509 //ofstream os("dbg.dat");
510 for (running_box->rewind();running_box->haveMore();
511 running_box->next()) {
512 const int ch=running_box->getChannel();
513 if (running_box->getNumberOfBoxPoints()>=minboxnchan)
514 processChannel(mask[ch] && (fabs(running_box->aboveMean()) >=
515 threshold*offline_variance), mask);
516 else processCurLine(mask); // just finish what was accumulated before
517
518 signs[ch]=getAboveMeanSign();
519 // os<<ch<<" "<<spectrum[ch]<<" "<<fabs(running_box->aboveMean())<<" "<<
520 // threshold*offline_variance<<endl;
521
522 const Float buf=running_box->aboveMean();
523 if (buf>0) signs[ch]=1;
524 else if (buf<0) signs[ch]=-1;
525 else if (buf==0) signs[ch]=0;
526 // os<<ch<<" "<<spectrum[ch]<<" "<<running_box->getLinMean()<<" "<<
527 // threshold*offline_variance<<endl;
528 }
529 if (lines.size())
530 searchForWings(lines,signs,mask,edge);
531 }
532 catch (const AipsError &ae) {
533 throw;
534 }
535 catch (const exception &ex) {
536 throw AipsError(String("LFAboveThreshold::findLines - STL error: ")+ex.what());
537 }
538}
539
540//
541///////////////////////////////////////////////////////////////////////////////
542
543///////////////////////////////////////////////////////////////////////////////
544//
545// LFLineListOperations::IntersectsWith - An auxiliary object function
546// to test whether two lines have a non-void intersection
547//
548
549
550// line1 - range of the first line: start channel and stop+1
551LFLineListOperations::IntersectsWith::IntersectsWith(const std::pair<int,int> &in_line1) :
552 line1(in_line1) {}
553
554
555// return true if line2 intersects with line1 with at least one
556// common channel, and false otherwise
557// line2 - range of the second line: start channel and stop+1
558bool LFLineListOperations::IntersectsWith::operator()(const std::pair<int,int> &line2)
559 const throw()
560{
561 if (line2.second<line1.first) return false; // line2 is at lower channels
562 if (line2.first>line1.second) return false; // line2 is at upper channels
563 return true; // line2 has an intersection or is adjacent to line1
564}
565
566//
567///////////////////////////////////////////////////////////////////////////////
568
569///////////////////////////////////////////////////////////////////////////////
570//
571// LFLineListOperations::BuildUnion - An auxiliary object function to build a union
572// of several lines to account for a possibility of merging the nearby lines
573//
574
575// set an initial line (can be a first line in the sequence)
576LFLineListOperations::BuildUnion::BuildUnion(const std::pair<int,int> &line1) :
577 temp_line(line1) {}
578
579// update temp_line with a union of temp_line and new_line
580// provided there is no gap between the lines
581void LFLineListOperations::BuildUnion::operator()(const std::pair<int,int> &new_line)
582 throw()
583{
584 if (new_line.first<temp_line.first) temp_line.first=new_line.first;
585 if (new_line.second>temp_line.second) temp_line.second=new_line.second;
586}
587
588// return the result (temp_line)
589const std::pair<int,int>& LFLineListOperations::BuildUnion::result() const throw()
590{
591 return temp_line;
592}
593
594//
595///////////////////////////////////////////////////////////////////////////////
596
597///////////////////////////////////////////////////////////////////////////////
598//
599// LFLineListOperations::LaterThan - An auxiliary object function to test whether a
600// specified line is at lower spectral channels (to preserve the order in
601// the line list)
602//
603
604// setup the line to compare with
605LFLineListOperations::LaterThan::LaterThan(const std::pair<int,int> &in_line1) :
606 line1(in_line1) {}
607
608// return true if line2 should be placed later than line1
609// in the ordered list (so, it is at greater channel numbers)
610bool LFLineListOperations::LaterThan::operator()(const std::pair<int,int> &line2)
611 const throw()
612{
613 if (line2.second<line1.first) return false; // line2 is at lower channels
614 if (line2.first>line1.second) return true; // line2 is at upper channels
615
616 // line2 intersects with line1. We should have no such situation in
617 // practice
618 return line2.first>line1.first;
619}
620
621//
622///////////////////////////////////////////////////////////////////////////////
623
624
625///////////////////////////////////////////////////////////////////////////////
626//
627// STLineFinder - a class for automated spectral line search
628//
629//
630
631STLineFinder::STLineFinder() throw() : edge(0,0)
632{
633 setOptions();
634}
635
636// set the parameters controlling algorithm
637// in_threshold a single channel threshold default is sqrt(3), which
638// means together with 3 minimum channels at least 3 sigma
639// detection criterion
640// For bad baseline shape, in_threshold may need to be
641// increased
642// in_min_nchan minimum number of channels above the threshold to report
643// a detection, default is 3
644// in_avg_limit perform the averaging of no more than in_avg_limit
645// adjacent channels to search for broad lines
646// Default is 8, but for a bad baseline shape this
647// parameter should be decreased (may be even down to a
648// minimum of 1 to disable this option) to avoid
649// confusing of baseline undulations with a real line.
650// Setting a very large value doesn't usually provide
651// valid detections.
652// in_box_size the box size for running mean calculation. Default is
653// 1./5. of the whole spectrum size
654void STLineFinder::setOptions(const casa::Float &in_threshold,
655 const casa::Int &in_min_nchan,
656 const casa::Int &in_avg_limit,
657 const casa::Float &in_box_size) throw()
658{
659 threshold=in_threshold;
660 min_nchan=in_min_nchan;
661 avg_limit=in_avg_limit;
662 box_size=in_box_size;
663}
664
665STLineFinder::~STLineFinder() throw(AipsError) {}
666
667// set scan to work with (in_scan parameter)
668void STLineFinder::setScan(const ScantableWrapper &in_scan) throw(AipsError)
669{
670 scan=in_scan.getCP();
671 AlwaysAssert(!scan.null(),AipsError);
672
673}
674
675// search for spectral lines. Number of lines found is returned
676// in_edge and in_mask control channel rejection for a given row
677// if in_edge has zero length, all channels chosen by mask will be used
678// if in_edge has one element only, it represents the number of
679// channels to drop from both sides of the spectrum
680// in_edge is introduced for convinience, although all functionality
681// can be achieved using a spectrum mask only
682int STLineFinder::findLines(const std::vector<bool> &in_mask,
683 const std::vector<int> &in_edge,
684 const casa::uInt &whichRow) throw(casa::AipsError)
685{
686 const int minboxnchan=4;
687 if (scan.null())
688 throw AipsError("STLineFinder::findLines - a scan should be set first,"
689 " use set_scan");
690
691 uInt nchan = scan->nchan(scan->getIF(whichRow));
692 // set up mask and edge rejection
693 // no mask given...
694 if (in_mask.size() == 0) {
695 mask = Vector<Bool>(nchan,True);
696 } else {
697 // use provided mask
698 mask=Vector<Bool>(in_mask);
699 }
700 if (mask.nelements()!=nchan)
701 throw AipsError("STLineFinder::findLines - in_scan and in_mask have different"
702 "number of spectral channels.");
703 cout << mask.nelements() << endl;
704 // number of elements in in_edge
705 if (in_edge.size()>2)
706 throw AipsError("STLineFinder::findLines - the length of the in_edge parameter"
707 "should not exceed 2");
708 if (!in_edge.size()) {
709 // all spectra, no rejection
710 edge.first=0;
711 edge.second=nchan;
712 } else {
713 edge.first=in_edge[0];
714 if (edge.first<0)
715 throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
716 "number of channels to drop");
717 if (edge.first>=nchan)
718 throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
719 if (in_edge.size()==2) {
720 edge.second=in_edge[1];
721 if (edge.second<0)
722 throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
723 "number of channels to drop");
724 edge.second=nchan-edge.second;
725 } else edge.second=nchan-edge.first;
726 if (edge.second<0 || (edge.first>=edge.second))
727 throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
728 }
729
730 //
731 int max_box_nchan=int(nchan*box_size); // number of channels in running
732 // box
733 if (max_box_nchan<2)
734 throw AipsError("STLineFinder::findLines - box_size is too small");
735
736 spectrum.resize();
737 spectrum = Vector<Float>(scan->getSpectrum(whichRow));
738
739 lines.resize(0); // search from the scratch
740 last_row_used=whichRow;
741 Vector<Bool> temp_mask(mask);
742
743 Bool first_pass=True;
744 Int avg_factor=1; // this number of adjacent channels is averaged together
745 // the total number of the channels is not altered
746 // instead, min_nchan is also scaled
747 // it helps to search for broad lines
748 Vector<Int> signs; // a buffer for signs of the value - mean quantity
749 // see LFAboveThreshold for details
750 // We need only signs resulted from last iteration
751 // because all previous values may be corrupted by the
752 // presence of spectral lines
753 while (true) {
754 // a buffer for new lines found at this iteration
755 std::list<pair<int,int> > new_lines;
756
757 try {
758 // line find algorithm
759 LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold);
760 lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan);
761 signs.resize(lfalg.getSigns().nelements());
762 signs=lfalg.getSigns();
763 first_pass=False;
764 if (!new_lines.size())
765 throw AipsError("spurious"); // nothing new - use the same
766 // code as for a real exception
767 }
768 catch(const AipsError &ae) {
769 if (first_pass) throw;
770 // nothing new - proceed to the next step of averaging, if any
771 // (to search for broad lines)
772 if (avg_factor>avg_limit) break; // averaging up to avg_limit
773 // adjacent channels,
774 // stop after that
775 avg_factor*=2; // twice as more averaging
776 subtractBaseline(temp_mask,9);
777 averageAdjacentChannels(temp_mask,avg_factor);
778 continue;
779 }
780 keepStrongestOnly(temp_mask,new_lines,max_box_nchan);
781 // update the list (lines) merging intervals, if necessary
782 addNewSearchResult(new_lines,lines);
783 // get a new mask
784 temp_mask=getMask();
785 }
786
787 // an additional search for wings because in the presence of very strong
788 // lines temporary mean used at each iteration will be higher than
789 // the true mean
790
791 if (lines.size())
792 LFLineListOperations::searchForWings(lines,signs,mask,edge);
793
794 return int(lines.size());
795}
796
797// auxiliary function to fit and subtract a polynomial from the current
798// spectrum. It uses the Fitter class. This action is required before
799// reducing the spectral resolution if the baseline shape is bad
800void STLineFinder::subtractBaseline(const casa::Vector<casa::Bool> &temp_mask,
801 const casa::Int &order) throw(casa::AipsError)
802{
803 AlwaysAssert(spectrum.nelements(),AipsError);
804 // use the fact that temp_mask excludes channels rejected at the edge
805 Fitter sdf;
806 std::vector<float> absc(spectrum.nelements());
807 for (Int i=0;i<absc.size();++i)
808 absc[i]=float(i)/float(spectrum.nelements());
809 std::vector<float> spec;
810 spectrum.tovector(spec);
811 std::vector<bool> std_mask;
812 temp_mask.tovector(std_mask);
813 sdf.setData(absc,spec,std_mask);
814 sdf.setExpression("poly",order);
815 if (!sdf.fit()) return; // fit failed, use old spectrum
816 spectrum=casa::Vector<casa::Float>(sdf.getResidual());
817}
818
819// auxiliary function to average adjacent channels and update the mask
820// if at least one channel involved in summation is masked, all
821// output channels will be masked. This function works with the
822// spectrum and edge fields of this class, but updates the mask
823// array specified, rather than the field of this class
824// boxsize - a number of adjacent channels to average
825void STLineFinder::averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update,
826 const casa::Int &boxsize)
827 throw(casa::AipsError)
828{
829 DebugAssert(mask2update.nelements()==spectrum.nelements(), AipsError);
830 DebugAssert(boxsize!=0,AipsError);
831
832 for (int n=edge.first;n<edge.second;n+=boxsize) {
833 DebugAssert(n<spectrum.nelements(),AipsError);
834 int nboxch=0; // number of channels currently in the box
835 Float mean=0; // buffer for mean calculations
836 for (int k=n;k<n+boxsize && k<edge.second;++k)
837 if (mask2update[k]) { // k is a valid channel
838 mean+=spectrum[k];
839 ++nboxch;
840 }
841 if (nboxch<boxsize) // mask these channels
842 for (int k=n;k<n+boxsize && k<edge.second;++k)
843 mask2update[k]=False;
844 else {
845 mean/=Float(boxsize);
846 for (int k=n;k<n+boxsize && k<edge.second;++k)
847 spectrum[k]=mean;
848 }
849 }
850}
851
852
853// get the mask to mask out all lines that have been found (default)
854// if invert=true, only channels belong to lines will be unmasked
855// Note: all channels originally masked by the input mask (in_mask
856// in setScan) or dropped out by the edge parameter (in_edge
857// in setScan) are still excluded regardless on the invert option
858std::vector<bool> STLineFinder::getMask(bool invert)
859 const throw(casa::AipsError)
860{
861 try {
862 if (scan.null())
863 throw AipsError("STLineFinder::getMask - a scan should be set first,"
864 " use set_scan followed by find_lines");
865 DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
866 /*
867 if (!lines.size())
868 throw AipsError("STLineFinder::getMask - one have to search for "
869 "lines first, use find_lines");
870 */
871 std::vector<bool> res_mask(mask.nelements());
872 // iterator through lines
873 std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
874 for (int ch=0;ch<res_mask.size();++ch)
875 if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
876 else if (!mask[ch]) res_mask[ch]=false;
877 else {
878 res_mask[ch]=!invert; // no line by default
879 if (cli==lines.end()) continue;
880 if (ch>=cli->first && ch<cli->second)
881 res_mask[ch]=invert; // this is a line
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// the sidelobes 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.