source: trunk/src/STLineFinder.cpp@ 2242

Last change on this file since 2242 was 2196, checked in by WataruKawasaki, 13 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): linefinder

Description: use LinearFit rather than NonlinearFitLM in STLineFinder.subtractBaseline(). this makes linefinder more than 10 times faster.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 52.6 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 2196 2011-06-21 10:14:05Z WataruKawasaki $
30//#---------------------------------------------------------------------------
31
32
33// ASAP
34#include "STLineFinder.h"
35#include "STFitter.h"
36#include "IndexedCompare.h"
37
38// STL
39#include <functional>
40#include <algorithm>
41#include <iostream>
42#include <fstream>
43
44using namespace asap;
45using namespace casa;
46using namespace std;
47
48namespace asap {
49
50///////////////////////////////////////////////////////////////////////////////
51//
52// RunningBox - a running box calculator. This class implements
53// iterations 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 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/median calculations.
114 // It adds a specified channel to the running 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 //
155 bool itsUseMedian; // true if median statistics is used
156 // to determine the noise level, otherwise
157 // it is the mean of the lowest 80% of deviations
158 // (default)
159 int itsNoiseSampleSize; // sample size used to estimate the noise statistics
160 // Negative value means the whole spectrum is used (default)
161public:
162
163 // set up the detection criterion
164 LFAboveThreshold(std::list<pair<int,int> > &in_lines,
165 int in_min_nchan = 3,
166 casa::Float in_threshold = 5,
167 bool use_median = false,
168 int noise_sample_size = -1) throw();
169 virtual ~LFAboveThreshold() throw();
170
171 // replace the detection criterion
172 void setCriterion(int in_min_nchan, casa::Float in_threshold) throw();
173
174 // return the array with signs of the value-current mean
175 // An element is +1 if value>mean, -1 if less, 0 if equal.
176 // This array is updated each time the findLines method is called and
177 // is used to search the line wings
178 const casa::Vector<Int>& getSigns() const throw();
179
180 // find spectral lines and add them into list
181 // if statholder is not NULL, the accumulate function of it will be
182 // called for each channel to save statistics
183 // spectrum, mask and edge - reference to the data
184 // max_box_nchan - number of channels in the running box
185 void findLines(const casa::Vector<casa::Float> &spectrum,
186 const casa::Vector<casa::Bool> &mask,
187 const std::pair<int,int> &edge,
188 int max_box_nchan) throw(casa::AipsError);
189
190protected:
191
192 // process a channel: update curline and is_detected before and
193 // add a new line to the list, if necessary using processCurLine()
194 // detect=true indicates that the current channel satisfies the criterion
195 void processChannel(Bool detect, const casa::Vector<casa::Bool> &mask)
196 throw(casa::AipsError);
197
198 // process the interval of channels stored in curline
199 // if it satisfies the criterion, add this interval as a new line
200 void processCurLine(const casa::Vector<casa::Bool> &mask)
201 throw(casa::AipsError);
202
203 // get the sign of runningBox->aboveMean(). The RunningBox pointer
204 // should be defined
205 casa::Int getAboveMeanSign() const throw();
206};
207
208//
209///////////////////////////////////////////////////////////////////////////////
210
211///////////////////////////////////////////////////////////////////////////////
212//
213// LFNoiseEstimator a helper class designed to estimate off-line variance
214// using statistics depending on the distribution of
215// values (e.g. like a median)
216//
217// Two statistics are supported: median and an average of
218// 80% of smallest values.
219//
220
221struct LFNoiseEstimator {
222 // construct an object
223 // size - maximum sample size. After a size number of elements is processed
224 // any new samples would cause the algorithm to drop the oldest samples in the
225 // buffer.
226 explicit LFNoiseEstimator(size_t size);
227
228 // add a new sample
229 // in - the new value
230 void add(float in);
231
232 // median of the distribution
233 float median() const;
234
235 // mean of lowest 80% of the samples
236 float meanLowest80Percent() const;
237
238 // return true if the buffer is full (i.e. statistics are representative)
239 inline bool filledToCapacity() const { return itsBufferFull;}
240
241protected:
242 // update cache of sorted indices
243 // (it is assumed that itsSampleNumber points to the newly
244 // replaced element)
245 void updateSortedCache() const;
246
247 // build sorted cache from the scratch
248 void buildSortedCache() const;
249
250 // number of samples accumulated so far
251 // (can be less than the buffer size)
252 size_t numberOfSamples() const;
253
254 // this helper method builds the cache if
255 // necessary using one of the methods
256 void fillCacheIfNecessary() const;
257
258private:
259 // buffer with samples (unsorted)
260 std::vector<float> itsVariances;
261 // current sample number (<=itsVariances.size())
262 size_t itsSampleNumber;
263 // true, if the buffer all values in the sample buffer are used
264 bool itsBufferFull;
265 // cached indices into vector of samples
266 mutable std::vector<size_t> itsSortedIndices;
267 // true if any of the statistics have been obtained at least
268 // once. This flag allows to implement a more efficient way of
269 // calculating statistics, if they are needed at once and not
270 // after each addition of a new element
271 mutable bool itsStatisticsAccessed;
272};
273
274//
275///////////////////////////////////////////////////////////////////////////////
276
277
278} // namespace asap
279
280///////////////////////////////////////////////////////////////////////////////
281//
282// LFNoiseEstimator a helper class designed to estimate off-line variance
283// using statistics depending on the distribution of
284// values (e.g. like a median)
285//
286// Two statistics are supported: median and an average of
287// 80% of smallest values.
288//
289
290// construct an object
291// size - maximum sample size. After a size number of elements is processed
292// any new samples would cause the algorithm to drop the oldest samples in the
293// buffer.
294LFNoiseEstimator::LFNoiseEstimator(size_t size) : itsVariances(size),
295 itsSampleNumber(0), itsBufferFull(false), itsSortedIndices(size),
296 itsStatisticsAccessed(false)
297{
298 AlwaysAssert(size>0,AipsError);
299}
300
301
302// add a new sample
303// in - the new value
304void LFNoiseEstimator::add(float in)
305{
306 if (isnan(in)) {
307 // normally it shouldn't happen
308 return;
309 }
310 itsVariances[itsSampleNumber] = in;
311
312 if (itsStatisticsAccessed) {
313 // only do element by element addition if on-the-fly
314 // statistics are needed
315 updateSortedCache();
316 }
317
318 // advance itsSampleNumber now
319 ++itsSampleNumber;
320 if (itsSampleNumber == itsVariances.size()) {
321 itsSampleNumber = 0;
322 itsBufferFull = true;
323 }
324 AlwaysAssert(itsSampleNumber<itsVariances.size(),AipsError);
325}
326
327// number of samples accumulated so far
328// (can be less than the buffer size)
329size_t LFNoiseEstimator::numberOfSamples() const
330{
331 // the number of samples accumulated so far may be less than the
332 // buffer size
333 const size_t nSamples = itsBufferFull ? itsVariances.size(): itsSampleNumber;
334 AlwaysAssert( (nSamples > 0) && (nSamples <= itsVariances.size()), AipsError);
335 return nSamples;
336}
337
338// this helper method builds the cache if
339// necessary using one of the methods
340void LFNoiseEstimator::fillCacheIfNecessary() const
341{
342 if (!itsStatisticsAccessed) {
343 if ((itsSampleNumber!=0) || itsBufferFull) {
344 // build the whole cache efficiently
345 buildSortedCache();
346 } else {
347 updateSortedCache();
348 }
349 itsStatisticsAccessed = true;
350 } // otherwise, it is updated in 'add' using on-the-fly method
351}
352
353// median of the distribution
354float LFNoiseEstimator::median() const
355{
356 fillCacheIfNecessary();
357 // the number of samples accumulated so far may be less than the
358 // buffer size
359 const size_t nSamples = numberOfSamples();
360 const size_t medSample = nSamples / 2;
361 AlwaysAssert(medSample < itsSortedIndices.size(), AipsError);
362 return itsVariances[itsSortedIndices[medSample]];
363}
364
365// mean of lowest 80% of the samples
366float LFNoiseEstimator::meanLowest80Percent() const
367{
368 fillCacheIfNecessary();
369 // the number of samples accumulated so far may be less than the
370 // buffer size
371 const size_t nSamples = numberOfSamples();
372 float result = 0;
373 size_t numpt=size_t(0.8*nSamples);
374 if (!numpt) {
375 numpt=nSamples; // no much else left,
376 // although it is very inaccurate
377 }
378 AlwaysAssert( (numpt > 0) && (numpt<itsSortedIndices.size()), AipsError);
379 for (size_t ch=0; ch<numpt; ++ch) {
380 result += itsVariances[itsSortedIndices[ch]];
381 }
382 result /= float(numpt);
383 return result;
384}
385
386// update cache of sorted indices
387// (it is assumed that itsSampleNumber points to the newly
388// replaced element)
389void LFNoiseEstimator::updateSortedCache() const
390{
391 // the number of samples accumulated so far may be less than the
392 // buffer size
393 const size_t nSamples = numberOfSamples();
394
395 if (itsBufferFull) {
396 // first find the index of the element which is being replaced
397 size_t index = nSamples;
398 for (size_t i=0; i<nSamples; ++i) {
399 AlwaysAssert(i < itsSortedIndices.size(), AipsError);
400 if (itsSortedIndices[i] == itsSampleNumber) {
401 index = i;
402 break;
403 }
404 }
405 AlwaysAssert( index < nSamples, AipsError);
406
407 const vector<size_t>::iterator indStart = itsSortedIndices.begin();
408 // merge this element with preceeding block first
409 if (index != 0) {
410 // merge indices on the basis of variances
411 inplace_merge(indStart,indStart+index,indStart+index+1,
412 indexedCompare<size_t>(itsVariances.begin()));
413 }
414 // merge with the following block
415 if (index + 1 != nSamples) {
416 // merge indices on the basis of variances
417 inplace_merge(indStart,indStart+index+1,indStart+nSamples,
418 indexedCompare<size_t>(itsVariances.begin()));
419 }
420 } else {
421 // itsSampleNumber is the index of the new element
422 AlwaysAssert(itsSampleNumber < itsSortedIndices.size(), AipsError);
423 itsSortedIndices[itsSampleNumber] = itsSampleNumber;
424 if (itsSampleNumber >= 1) {
425 // we have to place this new sample in
426 const vector<size_t>::iterator indStart = itsSortedIndices.begin();
427 // merge indices on the basis of variances
428 inplace_merge(indStart,indStart+itsSampleNumber,indStart+itsSampleNumber+1,
429 indexedCompare<size_t>(itsVariances.begin()));
430 }
431 }
432}
433
434// build sorted cache from the scratch
435void LFNoiseEstimator::buildSortedCache() const
436{
437 // the number of samples accumulated so far may be less than the
438 // buffer size
439 const size_t nSamples = numberOfSamples();
440 AlwaysAssert(nSamples <= itsSortedIndices.size(), AipsError);
441 for (size_t i=0; i<nSamples; ++i) {
442 itsSortedIndices[i]=i;
443 }
444
445 // sort indices, but check the array of variances
446 const vector<size_t>::iterator indStart = itsSortedIndices.begin();
447 stable_sort(indStart,indStart+nSamples, indexedCompare<size_t>(itsVariances.begin()));
448}
449
450//
451///////////////////////////////////////////////////////////////////////////////
452
453///////////////////////////////////////////////////////////////////////////////
454//
455// RunningBox - a running box calculator. This class implements
456// interations over the specified spectrum and calculates
457// running box filter statistics.
458//
459
460// set up the object with the references to actual data
461// and the number of channels in the running box
462RunningBox::RunningBox(const casa::Vector<casa::Float> &in_spectrum,
463 const casa::Vector<casa::Bool> &in_mask,
464 const std::pair<int,int> &in_edge,
465 int in_max_box_nchan) throw(AipsError) :
466 spectrum(in_spectrum), mask(in_mask), edge(in_edge),
467 max_box_nchan(in_max_box_nchan)
468{
469 rewind();
470}
471
472void RunningBox::rewind() throw(AipsError) {
473 // fill statistics for initial box
474 box_chan_cntr=0; // no channels are currently in the box
475 sumf=0.; // initialize statistics
476 sumf2=0.;
477 sumch=0.;
478 sumch2=0.;
479 sumfch=0.;
480 int initial_box_ch=edge.first;
481 for (;initial_box_ch<edge.second && box_chan_cntr<max_box_nchan;
482 ++initial_box_ch)
483 advanceRunningBox(initial_box_ch);
484
485 if (initial_box_ch==edge.second)
486 throw AipsError("RunningBox::rewind - too much channels are masked");
487
488 cur_channel=edge.first;
489 start_advance=initial_box_ch-max_box_nchan/2;
490}
491
492// access to the statistics
493const casa::Float& RunningBox::getLinMean() const throw(AipsError)
494{
495 DebugAssert(cur_channel<edge.second, AipsError);
496 if (need2recalculate) updateDerivativeStatistics();
497 return linmean;
498}
499
500const casa::Float& RunningBox::getLinVariance() const throw(AipsError)
501{
502 DebugAssert(cur_channel<edge.second, AipsError);
503 if (need2recalculate) updateDerivativeStatistics();
504 return linvariance;
505}
506
507casa::Float RunningBox::aboveMean() const throw(AipsError)
508{
509 DebugAssert(cur_channel<edge.second, AipsError);
510 if (need2recalculate) updateDerivativeStatistics();
511 return spectrum[cur_channel]-linmean;
512}
513
514int RunningBox::getChannel() const throw()
515{
516 return cur_channel;
517}
518
519// actual number of channels in the box (max_box_nchan, if no channels
520// are masked)
521int RunningBox::getNumberOfBoxPoints() const throw()
522{
523 return box_chan_cntr;
524}
525
526// supplementary function to control running mean/median calculations.
527// It adds a specified channel to the running box and
528// removes (ch-max_box_nchan+1)'th channel from there
529// Channels, for which the mask is false or index is beyond the
530// allowed range, are ignored
531void RunningBox::advanceRunningBox(int ch) throw(AipsError)
532{
533 if (ch>=edge.first && ch<edge.second)
534 if (mask[ch]) { // ch is a valid channel
535 ++box_chan_cntr;
536 sumf+=spectrum[ch];
537 sumf2+=square(spectrum[ch]);
538 sumch+=Float(ch);
539 sumch2+=square(Float(ch));
540 sumfch+=spectrum[ch]*Float(ch);
541 need2recalculate=True;
542 }
543 int ch2remove=ch-max_box_nchan;
544 if (ch2remove>=edge.first && ch2remove<edge.second)
545 if (mask[ch2remove]) { // ch2remove is a valid channel
546 --box_chan_cntr;
547 sumf-=spectrum[ch2remove];
548 sumf2-=square(spectrum[ch2remove]);
549 sumch-=Float(ch2remove);
550 sumch2-=square(Float(ch2remove));
551 sumfch-=spectrum[ch2remove]*Float(ch2remove);
552 need2recalculate=True;
553 }
554}
555
556// next channel
557void RunningBox::next() throw(AipsError)
558{
559 AlwaysAssert(cur_channel<edge.second,AipsError);
560 ++cur_channel;
561 if (cur_channel+max_box_nchan/2<edge.second && cur_channel>=start_advance)
562 advanceRunningBox(cur_channel+max_box_nchan/2); // update statistics
563}
564
565// checking whether there are still elements
566casa::Bool RunningBox::haveMore() const throw()
567{
568 return cur_channel<edge.second;
569}
570
571// calculate derivative statistics. This function is const, because
572// it updates the cache only
573void RunningBox::updateDerivativeStatistics() const throw(AipsError)
574{
575 AlwaysAssert(box_chan_cntr, AipsError);
576
577 Float mean=sumf/Float(box_chan_cntr);
578
579 // linear LSF formulae
580 Float meanch=sumch/Float(box_chan_cntr);
581 Float meanch2=sumch2/Float(box_chan_cntr);
582 if (meanch==meanch2 || box_chan_cntr<3) {
583 // vertical line in the spectrum, can't calculate linmean and linvariance
584 linmean=0.;
585 linvariance=0.;
586 } else {
587 Float coeff=(sumfch/Float(box_chan_cntr)-meanch*mean)/
588 (meanch2-square(meanch));
589 linmean=coeff*(Float(cur_channel)-meanch)+mean;
590 linvariance=sumf2/Float(box_chan_cntr)-square(mean)-
591 square(coeff)*(meanch2-square(meanch));
592 if (linvariance<0.) {
593 // this shouldn't happen normally, but could be due to round-off error
594 linvariance = 0;
595 } else {
596 linvariance = sqrt(linvariance);
597 }
598 }
599 need2recalculate=False;
600}
601
602
603//
604///////////////////////////////////////////////////////////////////////////////
605
606///////////////////////////////////////////////////////////////////////////////
607//
608// LFAboveThreshold - a running mean/median algorithm for line detection
609//
610//
611
612
613// set up the detection criterion
614LFAboveThreshold::LFAboveThreshold(std::list<pair<int,int> > &in_lines,
615 int in_min_nchan,
616 casa::Float in_threshold,
617 bool use_median,
618 int noise_sample_size) throw() :
619 min_nchan(in_min_nchan), threshold(in_threshold),
620 lines(in_lines), running_box(NULL), itsUseMedian(use_median),
621 itsNoiseSampleSize(noise_sample_size) {}
622
623LFAboveThreshold::~LFAboveThreshold() throw()
624{
625 if (running_box!=NULL) delete running_box;
626}
627
628// replace the detection criterion
629void LFAboveThreshold::setCriterion(int in_min_nchan, casa::Float in_threshold)
630 throw()
631{
632 min_nchan=in_min_nchan;
633 threshold=in_threshold;
634}
635
636// get the sign of runningBox->aboveMean(). The RunningBox pointer
637// should be defined
638casa::Int LFAboveThreshold::getAboveMeanSign() const throw()
639{
640 const Float buf=running_box->aboveMean();
641 if (buf>0) return 1;
642 if (buf<0) return -1;
643 return 0;
644}
645
646
647// process a channel: update cur_line and is_detected before and
648// add a new line to the list, if necessary
649void LFAboveThreshold::processChannel(Bool detect,
650 const casa::Vector<casa::Bool> &mask) throw(casa::AipsError)
651{
652 try {
653 if (is_detected_before) {
654 // we have to check that the current detection has the
655 // same sign of running_box->aboveMean
656 // otherwise it could be a spurious detection
657 if (last_sign && last_sign!=getAboveMeanSign())
658 detect=False;
659 }
660 if (detect) {
661 last_sign=getAboveMeanSign();
662 if (is_detected_before)
663 cur_line.second=running_box->getChannel()+1;
664 else {
665 is_detected_before=True;
666 cur_line.first=running_box->getChannel();
667 cur_line.second=running_box->getChannel()+1;
668 }
669 } else processCurLine(mask);
670 }
671 catch (const AipsError &ae) {
672 throw;
673 }
674 catch (const exception &ex) {
675 throw AipsError(String("LFAboveThreshold::processChannel - STL error: ")+ex.what());
676 }
677}
678
679// process the interval of channels stored in cur_line
680// if it satisfies the criterion, add this interval as a new line
681void LFAboveThreshold::processCurLine(const casa::Vector<casa::Bool> &mask)
682 throw(casa::AipsError)
683{
684 try {
685 if (is_detected_before) {
686 if (cur_line.second-cur_line.first>=min_nchan) {
687 // it was a detection. We need to change the list
688 Bool add_new_line=False;
689 if (lines.size()) {
690 for (int i=lines.back().second;i<cur_line.first;++i)
691 if (mask[i]) { // one valid channel in between
692 // means that we deal with a separate line
693 add_new_line=True;
694 break;
695 }
696 } else add_new_line=True;
697 if (add_new_line)
698 lines.push_back(cur_line);
699 else lines.back().second=cur_line.second;
700 }
701 is_detected_before=False;
702 }
703 }
704 catch (const AipsError &ae) {
705 throw;
706 }
707 catch (const exception &ex) {
708 throw AipsError(String("LFAboveThreshold::processCurLine - STL error: ")+ex.what());
709 }
710}
711
712// return the array with signs of the value-current mean
713// An element is +1 if value>mean, -1 if less, 0 if equal.
714// This array is updated each time the findLines method is called and
715// is used to search the line wings
716const casa::Vector<Int>& LFAboveThreshold::getSigns() const throw()
717{
718 return signs;
719}
720
721// find spectral lines and add them into list
722void LFAboveThreshold::findLines(const casa::Vector<casa::Float> &spectrum,
723 const casa::Vector<casa::Bool> &mask,
724 const std::pair<int,int> &edge,
725 int max_box_nchan)
726 throw(casa::AipsError)
727{
728 const int minboxnchan=4;
729 try {
730
731 if (running_box!=NULL) delete running_box;
732 running_box=new RunningBox(spectrum,mask,edge,max_box_nchan);
733
734 // determine the off-line variance first
735 // an assumption made: lines occupy a small part of the spectrum
736
737 const size_t noiseSampleSize = itsNoiseSampleSize<0 ? size_t(edge.second-edge.first) :
738 std::min(size_t(itsNoiseSampleSize), size_t(edge.second-edge.first));
739 DebugAssert(noiseSampleSize,AipsError);
740 const bool globalNoise = (size_t(edge.second - edge.first) == noiseSampleSize);
741 LFNoiseEstimator ne(noiseSampleSize);
742
743 for (;running_box->haveMore();running_box->next()) {
744 ne.add(running_box->getLinVariance());
745 if (ne.filledToCapacity()) {
746 break;
747 }
748 }
749
750 Float offline_variance = -1; // just a flag that it is unset
751
752 if (globalNoise) {
753 offline_variance = itsUseMedian ? ne.median() : ne.meanLowest80Percent();
754 }
755
756 // actual search algorithm
757 is_detected_before=False;
758
759 // initiate the signs array
760 signs.resize(spectrum.nelements());
761 signs=Vector<Int>(spectrum.nelements(),0);
762
763 //ofstream os("dbg.dat");
764 for (running_box->rewind();running_box->haveMore();
765 running_box->next()) {
766 const int ch=running_box->getChannel();
767 if (!globalNoise) {
768 // add a next point for a local noise estimate
769 ne.add(running_box->getLinVariance());
770 }
771 if (running_box->getNumberOfBoxPoints()>=minboxnchan) {
772 if (!globalNoise) {
773 offline_variance = itsUseMedian ? ne.median() : ne.meanLowest80Percent();
774 }
775 AlwaysAssert(offline_variance>0.,AipsError);
776 processChannel(mask[ch] && (fabs(running_box->aboveMean()) >=
777 threshold*offline_variance), mask);
778 } else processCurLine(mask); // just finish what was accumulated before
779
780 signs[ch]=getAboveMeanSign();
781 //os<<ch<<" "<<spectrum[ch]<<" "<<fabs(running_box->aboveMean())<<" "<<
782 //threshold*offline_variance<<endl;
783 }
784 if (lines.size())
785 searchForWings(lines,signs,mask,edge);
786 }
787 catch (const AipsError &ae) {
788 throw;
789 }
790 catch (const exception &ex) {
791 throw AipsError(String("LFAboveThreshold::findLines - STL error: ")+ex.what());
792 }
793}
794
795//
796///////////////////////////////////////////////////////////////////////////////
797
798///////////////////////////////////////////////////////////////////////////////
799//
800// LFLineListOperations::IntersectsWith - An auxiliary object function
801// to test whether two lines have a non-void intersection
802//
803
804
805// line1 - range of the first line: start channel and stop+1
806LFLineListOperations::IntersectsWith::IntersectsWith(const std::pair<int,int> &in_line1) :
807 line1(in_line1) {}
808
809
810// return true if line2 intersects with line1 with at least one
811// common channel, and false otherwise
812// line2 - range of the second line: start channel and stop+1
813bool LFLineListOperations::IntersectsWith::operator()(const std::pair<int,int> &line2)
814 const throw()
815{
816 if (line2.second<line1.first) return false; // line2 is at lower channels
817 if (line2.first>line1.second) return false; // line2 is at upper channels
818 return true; // line2 has an intersection or is adjacent to line1
819}
820
821//
822///////////////////////////////////////////////////////////////////////////////
823
824///////////////////////////////////////////////////////////////////////////////
825//
826// LFLineListOperations::BuildUnion - An auxiliary object function to build a union
827// of several lines to account for a possibility of merging the nearby lines
828//
829
830// set an initial line (can be a first line in the sequence)
831LFLineListOperations::BuildUnion::BuildUnion(const std::pair<int,int> &line1) :
832 temp_line(line1) {}
833
834// update temp_line with a union of temp_line and new_line
835// provided there is no gap between the lines
836void LFLineListOperations::BuildUnion::operator()(const std::pair<int,int> &new_line)
837 throw()
838{
839 if (new_line.first<temp_line.first) temp_line.first=new_line.first;
840 if (new_line.second>temp_line.second) temp_line.second=new_line.second;
841}
842
843// return the result (temp_line)
844const std::pair<int,int>& LFLineListOperations::BuildUnion::result() const throw()
845{
846 return temp_line;
847}
848
849//
850///////////////////////////////////////////////////////////////////////////////
851
852///////////////////////////////////////////////////////////////////////////////
853//
854// LFLineListOperations::LaterThan - An auxiliary object function to test whether a
855// specified line is at lower spectral channels (to preserve the order in
856// the line list)
857//
858
859// setup the line to compare with
860LFLineListOperations::LaterThan::LaterThan(const std::pair<int,int> &in_line1) :
861 line1(in_line1) {}
862
863// return true if line2 should be placed later than line1
864// in the ordered list (so, it is at greater channel numbers)
865bool LFLineListOperations::LaterThan::operator()(const std::pair<int,int> &line2)
866 const throw()
867{
868 if (line2.second<line1.first) return false; // line2 is at lower channels
869 if (line2.first>line1.second) return true; // line2 is at upper channels
870
871 // line2 intersects with line1. We should have no such situation in
872 // practice
873 return line2.first>line1.first;
874}
875
876//
877///////////////////////////////////////////////////////////////////////////////
878
879
880///////////////////////////////////////////////////////////////////////////////
881//
882// STLineFinder - a class for automated spectral line search
883//
884//
885
886STLineFinder::STLineFinder() throw() : edge(0,0)
887{
888 setOptions();
889}
890
891// set the parameters controlling algorithm
892// in_threshold a single channel threshold default is sqrt(3), which
893// means together with 3 minimum channels at least 3 sigma
894// detection criterion
895// For bad baseline shape, in_threshold may need to be
896// increased
897// in_min_nchan minimum number of channels above the threshold to report
898// a detection, default is 3
899// in_avg_limit perform the averaging of no more than in_avg_limit
900// adjacent channels to search for broad lines
901// Default is 8, but for a bad baseline shape this
902// parameter should be decreased (may be even down to a
903// minimum of 1 to disable this option) to avoid
904// confusing of baseline undulations with a real line.
905// Setting a very large value doesn't usually provide
906// valid detections.
907// in_box_size the box size for running mean/median calculation. Default is
908// 1./5. of the whole spectrum size
909// in_noise_box the box size for off-line noise estimation (if working with
910// local noise. Negative value means use global noise estimate
911// Default is -1 (i.e. estimate using the whole spectrum)
912// in_median true if median statistics is used as opposed to average of
913// the lowest 80% of deviations (default)
914void STLineFinder::setOptions(const casa::Float &in_threshold,
915 const casa::Int &in_min_nchan,
916 const casa::Int &in_avg_limit,
917 const casa::Float &in_box_size,
918 const casa::Float &in_noise_box,
919 const casa::Bool &in_median) throw()
920{
921 threshold=in_threshold;
922 min_nchan=in_min_nchan;
923 avg_limit=in_avg_limit;
924 box_size=in_box_size;
925 itsNoiseBox = in_noise_box;
926 itsUseMedian = in_median;
927
928 useScantable = true;
929}
930
931STLineFinder::~STLineFinder() throw(AipsError) {}
932
933// set scan to work with (in_scan parameter)
934void STLineFinder::setScan(const ScantableWrapper &in_scan) throw(AipsError)
935{
936 scan=in_scan.getCP();
937 AlwaysAssert(!scan.null(),AipsError);
938}
939
940// set spectrum data to work with. this is a method to allow linefinder work
941// without setting scantable for the purpose of using linefinder inside some
942// method in scantable class. (Dec 22, 2010 by W.Kawasaki)
943void STLineFinder::setData(const std::vector<float> &in_spectrum)
944{
945 spectrum = Vector<Float>(in_spectrum);
946 useScantable = false;
947}
948
949// search for spectral lines. Number of lines found is returned
950// in_edge and in_mask control channel rejection for a given row
951// if in_edge has zero length, all channels chosen by mask will be used
952// if in_edge has one element only, it represents the number of
953// channels to drop from both sides of the spectrum
954// in_edge is introduced for convinience, although all functionality
955// can be achieved using a spectrum mask only
956int STLineFinder::findLines(const std::vector<bool> &in_mask,
957 const std::vector<int> &in_edge,
958 const casa::uInt &whichRow) throw(casa::AipsError)
959{
960 if (useScantable && scan.null())
961 throw AipsError("STLineFinder::findLines - a scan should be set first,"
962 " use set_scan");
963
964 uInt nchan = useScantable ? scan->nchan(scan->getIF(whichRow)) : spectrum.nelements();
965 // set up mask and edge rejection
966 // no mask given...
967 if (in_mask.size() == 0) {
968 mask = Vector<Bool>(nchan,True);
969 } else {
970 // use provided mask
971 mask=Vector<Bool>(in_mask);
972 }
973 if (mask.nelements()!=nchan)
974 throw AipsError("STLineFinder::findLines - in_scan and in_mask, or in_spectrum "
975 "and in_mask have different number of spectral channels.");
976
977 // taking flagged channels into account
978 if (useScantable) {
979 vector<bool> flaggedChannels = scan->getMask(whichRow);
980 if (flaggedChannels.size()) {
981 // there is a mask set for this row
982 if (flaggedChannels.size() != mask.nelements()) {
983 throw AipsError("STLineFinder::findLines - internal inconsistency: number of "
984 "mask elements do not match the number of channels");
985 }
986 for (size_t ch = 0; ch<mask.nelements(); ++ch) {
987 mask[ch] &= flaggedChannels[ch];
988 }
989 }
990 }
991
992 // number of elements in in_edge
993 if (in_edge.size()>2)
994 throw AipsError("STLineFinder::findLines - the length of the in_edge parameter"
995 "should not exceed 2");
996 if (!in_edge.size()) {
997 // all spectra, no rejection
998 edge.first=0;
999 edge.second=nchan;
1000 } else {
1001 edge.first=in_edge[0];
1002 if (edge.first<0)
1003 throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
1004 "number of channels to drop");
1005 if (edge.first>=int(nchan))
1006 throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
1007 if (in_edge.size()==2) {
1008 edge.second=in_edge[1];
1009 if (edge.second<0)
1010 throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
1011 "number of channels to drop");
1012 edge.second=nchan-edge.second;
1013 } else edge.second=nchan-edge.first;
1014 if (edge.second<0 || (edge.first>=edge.second))
1015 throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
1016 }
1017
1018 //
1019 int max_box_nchan=int(nchan*box_size); // number of channels in running
1020 // box
1021 if (max_box_nchan<2)
1022 throw AipsError("STLineFinder::findLines - box_size is too small");
1023
1024 // number of elements in the sample for noise estimate
1025 const int noise_box = itsNoiseBox<0 ? -1 : int(nchan * itsNoiseBox);
1026
1027 if ((noise_box!= -1) and (noise_box<2))
1028 throw AipsError("STLineFinder::findLines - noise_box is supposed to be at least 2 elements");
1029
1030 if (useScantable) {
1031 spectrum.resize();
1032 spectrum = Vector<Float>(scan->getSpectrum(whichRow));
1033 }
1034
1035 lines.resize(0); // search from the scratch
1036 last_row_used=whichRow;
1037 Vector<Bool> temp_mask(mask);
1038
1039 Bool first_pass=True;
1040 Int avg_factor=1; // this number of adjacent channels is averaged together
1041 // the total number of the channels is not altered
1042 // instead, min_nchan is also scaled
1043 // it helps to search for broad lines
1044 Vector<Int> signs; // a buffer for signs of the value - mean quantity
1045 // see LFAboveThreshold for details
1046 // We need only signs resulted from last iteration
1047 // because all previous values may be corrupted by the
1048 // presence of spectral lines
1049 while (true) {
1050 // a buffer for new lines found at this iteration
1051 std::list<pair<int,int> > new_lines;
1052
1053 try {
1054 // line find algorithm
1055 LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold, itsUseMedian,noise_box);
1056 lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan);
1057 signs.resize(lfalg.getSigns().nelements());
1058 signs=lfalg.getSigns();
1059 first_pass=False;
1060 if (!new_lines.size())
1061 throw AipsError("spurious"); // nothing new - use the same
1062 // code as for a real exception
1063 }
1064 catch(const AipsError &ae) {
1065 if (first_pass) throw;
1066 // nothing new - proceed to the next step of averaging, if any
1067 // (to search for broad lines)
1068 if (avg_factor>=avg_limit) break; // averaging up to avg_limit
1069 // adjacent channels,
1070 // stop after that
1071 avg_factor*=2; // twice as more averaging
1072 subtractBaseline(temp_mask,9);
1073 averageAdjacentChannels(temp_mask,avg_factor);
1074 continue;
1075 }
1076 keepStrongestOnly(temp_mask,new_lines,max_box_nchan);
1077 // update the list (lines) merging intervals, if necessary
1078 addNewSearchResult(new_lines,lines);
1079 // get a new mask
1080 temp_mask=getMask();
1081 }
1082
1083 // an additional search for wings because in the presence of very strong
1084 // lines temporary mean used at each iteration will be higher than
1085 // the true mean
1086
1087 if (lines.size())
1088 LFLineListOperations::searchForWings(lines,signs,mask,edge);
1089
1090 return int(lines.size());
1091}
1092
1093// auxiliary function to fit and subtract a polynomial from the current
1094// spectrum. It uses the Fitter class. This action is required before
1095// reducing the spectral resolution if the baseline shape is bad
1096void STLineFinder::subtractBaseline(const casa::Vector<casa::Bool> &temp_mask,
1097 const casa::Int &order) throw(casa::AipsError)
1098{
1099 AlwaysAssert(spectrum.nelements(),AipsError);
1100 // use the fact that temp_mask excludes channels rejected at the edge
1101 Fitter sdf;
1102 std::vector<float> absc(spectrum.nelements());
1103 for (unsigned int i=0;i<absc.size();++i)
1104 absc[i]=float(i)/float(spectrum.nelements());
1105 std::vector<float> spec;
1106 spectrum.tovector(spec);
1107 std::vector<bool> std_mask;
1108 temp_mask.tovector(std_mask);
1109 sdf.setData(absc,spec,std_mask);
1110 sdf.setExpression("poly",order);
1111 if (!sdf.lfit()) return; // fit failed, use old spectrum
1112 spectrum=casa::Vector<casa::Float>(sdf.getResidual());
1113}
1114
1115// auxiliary function to average adjacent channels and update the mask
1116// if at least one channel involved in summation is masked, all
1117// output channels will be masked. This function works with the
1118// spectrum and edge fields of this class, but updates the mask
1119// array specified, rather than the field of this class
1120// boxsize - a number of adjacent channels to average
1121void STLineFinder::averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update,
1122 const casa::Int &boxsize)
1123 throw(casa::AipsError)
1124{
1125 DebugAssert(mask2update.nelements()==spectrum.nelements(), AipsError);
1126 DebugAssert(boxsize!=0,AipsError);
1127
1128 for (int n=edge.first;n<edge.second;n+=boxsize) {
1129 DebugAssert(n<spectrum.nelements(),AipsError);
1130 int nboxch=0; // number of channels currently in the box
1131 Float mean=0; // buffer for mean calculations
1132 for (int k=n;k<n+boxsize && k<edge.second;++k)
1133 if (mask2update[k]) { // k is a valid channel
1134 mean+=spectrum[k];
1135 ++nboxch;
1136 }
1137 if (nboxch<boxsize) // mask these channels
1138 for (int k=n;k<n+boxsize && k<edge.second;++k)
1139 mask2update[k]=False;
1140 else {
1141 mean/=Float(boxsize);
1142 for (int k=n;k<n+boxsize && k<edge.second;++k)
1143 spectrum[k]=mean;
1144 }
1145 }
1146}
1147
1148
1149// get the mask to mask out all lines that have been found (default)
1150// if invert=true, only channels belong to lines will be unmasked
1151// Note: all channels originally masked by the input mask (in_mask
1152// in setScan) or dropped out by the edge parameter (in_edge
1153// in setScan) are still excluded regardless on the invert option
1154std::vector<bool> STLineFinder::getMask(bool invert)
1155 const throw(casa::AipsError)
1156{
1157 try {
1158 if (useScantable) {
1159 if (scan.null())
1160 throw AipsError("STLineFinder::getMask - a scan should be set first,"
1161 " use set_scan followed by find_lines");
1162 DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
1163 }
1164 /*
1165 if (!lines.size())
1166 throw AipsError("STLineFinder::getMask - one have to search for "
1167 "lines first, use find_lines");
1168 */
1169 std::vector<bool> res_mask(mask.nelements());
1170 // iterator through lines
1171 std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
1172 for (int ch=0;ch<int(res_mask.size());++ch) {
1173 if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
1174 else if (!mask[ch]) res_mask[ch]=false;
1175 else {
1176 res_mask[ch]=!invert; // no line by default
1177 if (cli!=lines.end())
1178 if (ch>=cli->first && ch<cli->second)
1179 res_mask[ch]=invert; // this is a line
1180 }
1181 if (cli!=lines.end())
1182 if (ch>=cli->second)
1183 ++cli; // next line in the list
1184 }
1185 return res_mask;
1186 }
1187 catch (const AipsError &ae) {
1188 throw;
1189 }
1190 catch (const exception &ex) {
1191 throw AipsError(String("STLineFinder::getMask - STL error: ")+ex.what());
1192 }
1193}
1194
1195// get range for all lines found. The same units as used in the scan
1196// will be returned (e.g. velocity instead of channels).
1197std::vector<double> STLineFinder::getLineRanges()
1198 const throw(casa::AipsError)
1199{
1200 std::vector<double> vel;
1201 if (useScantable) {
1202 // convert to required abscissa units
1203 vel = scan->getAbcissa(last_row_used);
1204 } else {
1205 for (uInt i = 0; i < spectrum.nelements(); ++i)
1206 vel.push_back((double)i);
1207 }
1208 std::vector<int> ranges=getLineRangesInChannels();
1209 std::vector<double> res(ranges.size());
1210
1211 std::vector<int>::const_iterator cri=ranges.begin();
1212 std::vector<double>::iterator outi=res.begin();
1213 for (;cri!=ranges.end() && outi!=res.end();++cri,++outi)
1214 if (uInt(*cri)>=vel.size())
1215 throw AipsError("STLineFinder::getLineRanges - getAbcissa provided less channels than reqired");
1216 else *outi=vel[*cri];
1217 return res;
1218}
1219
1220// The same as getLineRanges, but channels are always used to specify
1221// the range
1222std::vector<int> STLineFinder::getLineRangesInChannels()
1223 const throw(casa::AipsError)
1224{
1225 try {
1226 if (useScantable) {
1227 if (scan.null())
1228 throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first,"
1229 " use set_scan followed by find_lines");
1230 DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
1231 }
1232
1233 if (!lines.size())
1234 throw AipsError("STLineFinder::getLineRangesInChannels - one have to search for "
1235 "lines first, use find_lines");
1236
1237 std::vector<int> res(2*lines.size());
1238 // iterator through lines & result
1239 std::list<std::pair<int,int> >::const_iterator cli = lines.begin();
1240 std::vector<int>::iterator ri = res.begin();
1241 for (; cli != lines.end() && ri != res.end(); ++cli,++ri) {
1242 *ri = cli->first;
1243 if (++ri != res.end())
1244 *ri = cli->second - 1;
1245 }
1246 return res;
1247 } catch (const AipsError &ae) {
1248 throw;
1249 } catch (const exception &ex) {
1250 throw AipsError(String("STLineFinder::getLineRanges - STL error: ") + ex.what());
1251 }
1252}
1253
1254
1255
1256// an auxiliary function to remove all lines from the list, except the
1257// strongest one (by absolute value). If the lines removed are real,
1258// they will be find again at the next iteration. This approach
1259// increases the number of iterations required, but is able to remove
1260// spurious detections likely to occur near strong lines.
1261// Later a better criterion may be implemented, e.g.
1262// taking into consideration the brightness of different lines. Now
1263// use the simplest solution
1264// temp_mask - mask to work with (may be different from original mask as
1265// the lines previously found may be masked)
1266// lines2update - a list of lines to work with
1267// nothing will be done if it is empty
1268// max_box_nchan - channels in the running box for baseline filtering
1269void STLineFinder::keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask,
1270 std::list<std::pair<int, int> > &lines2update,
1271 int max_box_nchan)
1272 throw (casa::AipsError)
1273{
1274 try {
1275 if (!lines2update.size()) return; // ignore an empty list
1276
1277 // current line
1278 std::list<std::pair<int,int> >::iterator li=lines2update.begin();
1279 // strongest line
1280 std::list<std::pair<int,int> >::iterator strongli=lines2update.begin();
1281 // the flux (absolute value) of the strongest line
1282 Float peak_flux=-1; // negative value - a flag showing uninitialized
1283 // value
1284 // the algorithm below relies on the list being ordered
1285 Float tmp_flux=-1; // a temporary peak
1286 for (RunningBox running_box(spectrum,temp_mask,edge,max_box_nchan);
1287 running_box.haveMore(); running_box.next()) {
1288
1289 if (li==lines2update.end()) break; // no more lines
1290 const int ch=running_box.getChannel();
1291 if (ch>=li->first && ch<li->second)
1292 if (temp_mask[ch] && tmp_flux<fabs(running_box.aboveMean()))
1293 tmp_flux=fabs(running_box.aboveMean());
1294 if (ch==li->second-1) {
1295 if (peak_flux<tmp_flux) { // if peak_flux=-1, this condition
1296 peak_flux=tmp_flux; // will be satisfied
1297 strongli=li;
1298 }
1299 ++li;
1300 tmp_flux=-1;
1301 }
1302 }
1303 std::list<std::pair<int,int> > res;
1304 res.splice(res.end(),lines2update,strongli);
1305 lines2update.clear();
1306 lines2update.splice(lines2update.end(),res);
1307 }
1308 catch (const AipsError &ae) {
1309 throw;
1310 }
1311 catch (const exception &ex) {
1312 throw AipsError(String("STLineFinder::keepStrongestOnly - STL error: ")+ex.what());
1313 }
1314
1315}
1316
1317//
1318///////////////////////////////////////////////////////////////////////////////
1319
1320
1321///////////////////////////////////////////////////////////////////////////////
1322//
1323// LFLineListOperations - a class incapsulating operations with line lists
1324// The LF prefix stands for Line Finder
1325//
1326
1327// concatenate two lists preserving the order. If two lines appear to
1328// be adjacent, they are joined into the new one
1329void LFLineListOperations::addNewSearchResult(const std::list<pair<int, int> > &newlines,
1330 std::list<std::pair<int, int> > &lines_list)
1331 throw(AipsError)
1332{
1333 try {
1334 for (std::list<pair<int,int> >::const_iterator cli=newlines.begin();
1335 cli!=newlines.end();++cli) {
1336
1337 // the first item, which has a non-void intersection or touches
1338 // the new line
1339 std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(),
1340 lines_list.end(), IntersectsWith(*cli));
1341 // the last such item
1342 std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg,
1343 lines_list.end(), not1(IntersectsWith(*cli)));
1344
1345 // extract all lines which intersect or touch a new one into
1346 // a temporary buffer. This may invalidate the iterators
1347 // line_buffer may be empty, if no lines intersects with a new
1348 // one.
1349 std::list<pair<int,int> > lines_buffer;
1350 lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end);
1351
1352 // build a union of all intersecting lines
1353 pair<int,int> union_line=for_each(lines_buffer.begin(),
1354 lines_buffer.end(),BuildUnion(*cli)).result();
1355
1356 // search for a right place for the new line (union_line) and add
1357 std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(),
1358 lines_list.end(), LaterThan(union_line));
1359 lines_list.insert(pos2insert,union_line);
1360 }
1361 }
1362 catch (const AipsError &ae) {
1363 throw;
1364 }
1365 catch (const exception &ex) {
1366 throw AipsError(String("LFLineListOperations::addNewSearchResult - STL error: ")+ex.what());
1367 }
1368}
1369
1370// extend all line ranges to the point where a value stored in the
1371// specified vector changes (e.g. value-mean change its sign)
1372// This operation is necessary to include line wings, which are below
1373// the detection threshold. If lines becomes adjacent, they are
1374// merged together. Any masked channel stops the extension
1375void LFLineListOperations::searchForWings(std::list<std::pair<int, int> > &newlines,
1376 const casa::Vector<casa::Int> &signs,
1377 const casa::Vector<casa::Bool> &mask,
1378 const std::pair<int,int> &edge) throw(casa::AipsError)
1379{
1380 try {
1381 for (std::list<pair<int,int> >::iterator li=newlines.begin();
1382 li!=newlines.end();++li) {
1383 // update the left hand side
1384 for (int n=li->first-1;n>=edge.first;--n) {
1385 if (!mask[n]) break;
1386 if (signs[n]==signs[li->first] && signs[li->first])
1387 li->first=n;
1388 else break;
1389 }
1390 // update the right hand side
1391 for (int n=li->second;n<edge.second;++n) {
1392 if (!mask[n]) break;
1393 if (signs[n]==signs[li->second-1] && signs[li->second-1])
1394 li->second=n;
1395 else break;
1396 }
1397 }
1398 // need to search for possible mergers.
1399 std::list<std::pair<int, int> > result_buffer;
1400 addNewSearchResult(newlines,result_buffer);
1401 newlines.clear();
1402 newlines.splice(newlines.end(),result_buffer);
1403 }
1404 catch (const AipsError &ae) {
1405 throw;
1406 }
1407 catch (const exception &ex) {
1408 throw AipsError(String("LFLineListOperations::extendLines - STL error: ")+ex.what());
1409 }
1410}
1411
1412//
1413///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.