source: branches/Release-1-fixes/src/SDLineFinder.cc@ 1418

Last change on this file since 1418 was 551, checked in by vor010, 20 years ago

Small improvements in the line finder. Most probably users will not notice them at all. Line wings are now detected better
for strong lines.

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