source: branches/hpc34/src/STLineFinder.cpp@ 2790

Last change on this file since 2790 was 2580, checked in by ShinnosukeKawakami, 12 years ago

hpc33 merged asap-trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 52.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 2580 2012-06-28 04:22:10Z ShinnosukeKawakami $
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), err("spurious")
887{
888 useScantable = true;
889 setOptions();
890}
891
892// set the parameters controlling algorithm
893// in_threshold a single channel threshold default is sqrt(3), which
894// means together with 3 minimum channels at least 3 sigma
895// detection criterion
896// For bad baseline shape, in_threshold may need to be
897// increased
898// in_min_nchan minimum number of channels above the threshold to report
899// a detection, default is 3
900// in_avg_limit perform the averaging of no more than in_avg_limit
901// adjacent channels to search for broad lines
902// Default is 8, but for a bad baseline shape this
903// parameter should be decreased (may be even down to a
904// minimum of 1 to disable this option) to avoid
905// confusing of baseline undulations with a real line.
906// Setting a very large value doesn't usually provide
907// valid detections.
908// in_box_size the box size for running mean/median calculation. Default is
909// 1./5. of the whole spectrum size
910// in_noise_box the box size for off-line noise estimation (if working with
911// local noise. Negative value means use global noise estimate
912// Default is -1 (i.e. estimate using the whole spectrum)
913// in_median true if median statistics is used as opposed to average of
914// the lowest 80% of deviations (default)
915void STLineFinder::setOptions(const casa::Float &in_threshold,
916 const casa::Int &in_min_nchan,
917 const casa::Int &in_avg_limit,
918 const casa::Float &in_box_size,
919 const casa::Float &in_noise_box,
920 const casa::Bool &in_median) throw()
921{
922 threshold=in_threshold;
923 min_nchan=in_min_nchan;
924 avg_limit=in_avg_limit;
925 box_size=in_box_size;
926 itsNoiseBox = in_noise_box;
927 itsUseMedian = in_median;
928}
929
930STLineFinder::~STLineFinder() throw(AipsError) {}
931
932// set scan to work with (in_scan parameter)
933void STLineFinder::setScan(const ScantableWrapper &in_scan) throw(AipsError)
934{
935 scan=in_scan.getCP();
936 AlwaysAssert(!scan.null(),AipsError);
937 useScantable = true;
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 spectrum.assign( Vector<Float>(in_spectrum) );
947 useScantable = false;
948}
949
950// search for spectral lines. Number of lines found is returned
951// in_edge and in_mask control channel rejection for a given row
952// if in_edge has zero length, all channels chosen by mask will be used
953// if in_edge has one element only, it represents the number of
954// channels to drop from both sides of the spectrum
955// in_edge is introduced for convinience, although all functionality
956// can be achieved using a spectrum mask only
957int STLineFinder::findLines(const std::vector<bool> &in_mask,
958 const std::vector<int> &in_edge,
959 const casa::uInt &whichRow) throw(casa::AipsError)
960{
961 if (useScantable && scan.null())
962 throw AipsError("STLineFinder::findLines - a scan should be set first,"
963 " use set_scan");
964
965 uInt nchan = useScantable ? scan->nchan(scan->getIF(whichRow)) : spectrum.nelements();
966 // set up mask and edge rejection
967 // no mask given...
968 if (in_mask.size() == 0) {
969 //mask = Vector<Bool>(nchan,True);
970 mask.assign( Vector<Bool>(nchan,True) );
971 } else {
972 // use provided mask
973 //mask=Vector<Bool>(in_mask);
974 mask.assign( Vector<Bool>(in_mask) );
975 }
976 if (mask.nelements()!=nchan)
977 throw AipsError("STLineFinder::findLines - in_scan and in_mask, or in_spectrum "
978 "and in_mask have different number of spectral channels.");
979
980 // taking flagged channels into account
981 if (useScantable) {
982 vector<bool> flaggedChannels = scan->getMask(whichRow);
983 if (flaggedChannels.size()) {
984 // there is a mask set for this row
985 if (flaggedChannels.size() != mask.nelements()) {
986 throw AipsError("STLineFinder::findLines - internal inconsistency: number of "
987 "mask elements do not match the number of channels");
988 }
989 for (size_t ch = 0; ch<mask.nelements(); ++ch) {
990 mask[ch] &= flaggedChannels[ch];
991 }
992 }
993 }
994
995 // number of elements in in_edge
996 if (in_edge.size()>2)
997 throw AipsError("STLineFinder::findLines - the length of the in_edge parameter"
998 "should not exceed 2");
999 if (!in_edge.size()) {
1000 // all spectra, no rejection
1001 edge.first=0;
1002 edge.second=nchan;
1003 } else {
1004 edge.first=in_edge[0];
1005 if (edge.first<0)
1006 throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
1007 "number of channels to drop");
1008 if (edge.first>=int(nchan))
1009 throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
1010 if (in_edge.size()==2) {
1011 edge.second=in_edge[1];
1012 if (edge.second<0)
1013 throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
1014 "number of channels to drop");
1015 edge.second=nchan-edge.second;
1016 } else edge.second=nchan-edge.first;
1017 if (edge.second<0 || (edge.first>=edge.second))
1018 throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
1019 }
1020
1021 //
1022 int max_box_nchan=int(nchan*box_size); // number of channels in running
1023 // box
1024 if (max_box_nchan<2)
1025 throw AipsError("STLineFinder::findLines - box_size is too small");
1026
1027 // number of elements in the sample for noise estimate
1028 const int noise_box = itsNoiseBox<0 ? -1 : int(nchan * itsNoiseBox);
1029
1030 if ((noise_box!= -1) and (noise_box<2))
1031 throw AipsError("STLineFinder::findLines - noise_box is supposed to be at least 2 elements");
1032
1033 if (useScantable) {
1034 spectrum.resize();
1035 spectrum = Vector<Float>(scan->getSpectrum(whichRow));
1036 }
1037
1038 lines.resize(0); // search from the scratch
1039 last_row_used=whichRow;
1040 Vector<Bool> temp_mask(mask);
1041
1042 Bool first_pass=True;
1043 Int avg_factor=1; // this number of adjacent channels is averaged together
1044 // the total number of the channels is not altered
1045 // instead, min_nchan is also scaled
1046 // it helps to search for broad lines
1047 Vector<Int> signs; // a buffer for signs of the value - mean quantity
1048 // see LFAboveThreshold for details
1049 // We need only signs resulted from last iteration
1050 // because all previous values may be corrupted by the
1051 // presence of spectral lines
1052
1053 while (true) {
1054 // a buffer for new lines found at this iteration
1055 std::list<pair<int,int> > new_lines;
1056
1057 try {
1058 // line find algorithm
1059 LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold, itsUseMedian,noise_box);
1060 lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan);
1061 signs.resize(lfalg.getSigns().nelements());
1062 signs=lfalg.getSigns();
1063 first_pass=False;
1064 if (!new_lines.size())
1065// throw AipsError("spurious"); // nothing new - use the same
1066// // code as for a real exception
1067 throw err; // nothing new - use the same
1068 // code as for a real exception
1069 }
1070 catch(const AipsError &ae) {
1071 if (first_pass) throw;
1072 // nothing new - proceed to the next step of averaging, if any
1073 // (to search for broad lines)
1074 if (avg_factor>=avg_limit) break; // averaging up to avg_limit
1075 // adjacent channels,
1076 // stop after that
1077 avg_factor*=2; // twice as more averaging
1078 subtractBaseline(temp_mask,9);
1079 averageAdjacentChannels(temp_mask,avg_factor);
1080 continue;
1081 }
1082 keepStrongestOnly(temp_mask,new_lines,max_box_nchan);
1083 // update the list (lines) merging intervals, if necessary
1084 addNewSearchResult(new_lines,lines);
1085 // get a new mask
1086 temp_mask=getMask();
1087 }
1088
1089 // an additional search for wings because in the presence of very strong
1090 // lines temporary mean used at each iteration will be higher than
1091 // the true mean
1092
1093 if (lines.size())
1094 LFLineListOperations::searchForWings(lines,signs,mask,edge);
1095
1096 return int(lines.size());
1097}
1098
1099// auxiliary function to fit and subtract a polynomial from the current
1100// spectrum. It uses the Fitter class. This action is required before
1101// reducing the spectral resolution if the baseline shape is bad
1102void STLineFinder::subtractBaseline(const casa::Vector<casa::Bool> &temp_mask,
1103 const casa::Int &order) throw(casa::AipsError)
1104{
1105 AlwaysAssert(spectrum.nelements(),AipsError);
1106 // use the fact that temp_mask excludes channels rejected at the edge
1107 Fitter sdf;
1108 std::vector<float> absc(spectrum.nelements());
1109 for (unsigned int i=0;i<absc.size();++i)
1110 absc[i]=float(i)/float(spectrum.nelements());
1111 std::vector<float> spec;
1112 spectrum.tovector(spec);
1113 std::vector<bool> std_mask;
1114 temp_mask.tovector(std_mask);
1115 sdf.setData(absc,spec,std_mask);
1116 sdf.setExpression("poly",order);
1117 if (!sdf.lfit()) return; // fit failed, use old spectrum
1118 spectrum=casa::Vector<casa::Float>(sdf.getResidual());
1119}
1120
1121// auxiliary function to average adjacent channels and update the mask
1122// if at least one channel involved in summation is masked, all
1123// output channels will be masked. This function works with the
1124// spectrum and edge fields of this class, but updates the mask
1125// array specified, rather than the field of this class
1126// boxsize - a number of adjacent channels to average
1127void STLineFinder::averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update,
1128 const casa::Int &boxsize)
1129 throw(casa::AipsError)
1130{
1131 DebugAssert(mask2update.nelements()==spectrum.nelements(), AipsError);
1132 DebugAssert(boxsize!=0,AipsError);
1133
1134 for (int n=edge.first;n<edge.second;n+=boxsize) {
1135 DebugAssert(n<spectrum.nelements(),AipsError);
1136 int nboxch=0; // number of channels currently in the box
1137 Float mean=0; // buffer for mean calculations
1138 for (int k=n;k<n+boxsize && k<edge.second;++k)
1139 if (mask2update[k]) { // k is a valid channel
1140 mean+=spectrum[k];
1141 ++nboxch;
1142 }
1143 if (nboxch<boxsize) // mask these channels
1144 for (int k=n;k<n+boxsize && k<edge.second;++k)
1145 mask2update[k]=False;
1146 else {
1147 mean/=Float(boxsize);
1148 for (int k=n;k<n+boxsize && k<edge.second;++k)
1149 spectrum[k]=mean;
1150 }
1151 }
1152}
1153
1154
1155// get the mask to mask out all lines that have been found (default)
1156// if invert=true, only channels belong to lines will be unmasked
1157// Note: all channels originally masked by the input mask (in_mask
1158// in setScan) or dropped out by the edge parameter (in_edge
1159// in setScan) are still excluded regardless on the invert option
1160std::vector<bool> STLineFinder::getMask(bool invert)
1161 const throw(casa::AipsError)
1162{
1163 try {
1164 if (useScantable) {
1165 if (scan.null())
1166 throw AipsError("STLineFinder::getMask - a scan should be set first,"
1167 " use set_scan followed by find_lines");
1168 DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
1169 }
1170 /*
1171 if (!lines.size())
1172 throw AipsError("STLineFinder::getMask - one have to search for "
1173 "lines first, use find_lines");
1174 */
1175 std::vector<bool> res_mask(mask.nelements());
1176 // iterator through lines
1177 std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
1178 for (int ch=0;ch<int(res_mask.size());++ch) {
1179 if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
1180 else if (!mask[ch]) res_mask[ch]=false;
1181 else {
1182 res_mask[ch]=!invert; // no line by default
1183 if (cli!=lines.end())
1184 if (ch>=cli->first && ch<cli->second)
1185 res_mask[ch]=invert; // this is a line
1186 }
1187 if (cli!=lines.end())
1188 if (ch>=cli->second)
1189 ++cli; // next line in the list
1190 }
1191 return res_mask;
1192 }
1193 catch (const AipsError &ae) {
1194 throw;
1195 }
1196 catch (const exception &ex) {
1197 throw AipsError(String("STLineFinder::getMask - STL error: ")+ex.what());
1198 }
1199}
1200
1201// get range for all lines found. The same units as used in the scan
1202// will be returned (e.g. velocity instead of channels).
1203std::vector<double> STLineFinder::getLineRanges()
1204 const throw(casa::AipsError)
1205{
1206 std::vector<double> vel;
1207 if (useScantable) {
1208 // convert to required abscissa units
1209 vel = scan->getAbcissa(last_row_used);
1210 } else {
1211 for (uInt i = 0; i < spectrum.nelements(); ++i)
1212 vel.push_back((double)i);
1213 }
1214 std::vector<int> ranges=getLineRangesInChannels();
1215 std::vector<double> res(ranges.size());
1216
1217 std::vector<int>::const_iterator cri=ranges.begin();
1218 std::vector<double>::iterator outi=res.begin();
1219 for (;cri!=ranges.end() && outi!=res.end();++cri,++outi)
1220 if (uInt(*cri)>=vel.size())
1221 throw AipsError("STLineFinder::getLineRanges - getAbcissa provided less channels than reqired");
1222 else *outi=vel[*cri];
1223 return res;
1224}
1225
1226// The same as getLineRanges, but channels are always used to specify
1227// the range
1228std::vector<int> STLineFinder::getLineRangesInChannels()
1229 const throw(casa::AipsError)
1230{
1231 try {
1232 if (useScantable) {
1233 if (scan.null())
1234 throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first,"
1235 " use set_scan followed by find_lines");
1236 DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
1237 }
1238
1239 if (!lines.size())
1240 throw AipsError("STLineFinder::getLineRangesInChannels - one have to search for "
1241 "lines first, use find_lines");
1242
1243 std::vector<int> res(2*lines.size());
1244 // iterator through lines & result
1245 std::list<std::pair<int,int> >::const_iterator cli = lines.begin();
1246 std::vector<int>::iterator ri = res.begin();
1247 for (; cli != lines.end() && ri != res.end(); ++cli,++ri) {
1248 *ri = cli->first;
1249 if (++ri != res.end())
1250 *ri = cli->second - 1;
1251 }
1252 return res;
1253 } catch (const AipsError &ae) {
1254 throw;
1255 } catch (const exception &ex) {
1256 throw AipsError(String("STLineFinder::getLineRanges - STL error: ") + ex.what());
1257 }
1258}
1259
1260
1261
1262// an auxiliary function to remove all lines from the list, except the
1263// strongest one (by absolute value). If the lines removed are real,
1264// they will be find again at the next iteration. This approach
1265// increases the number of iterations required, but is able to remove
1266// spurious detections likely to occur near strong lines.
1267// Later a better criterion may be implemented, e.g.
1268// taking into consideration the brightness of different lines. Now
1269// use the simplest solution
1270// temp_mask - mask to work with (may be different from original mask as
1271// the lines previously found may be masked)
1272// lines2update - a list of lines to work with
1273// nothing will be done if it is empty
1274// max_box_nchan - channels in the running box for baseline filtering
1275void STLineFinder::keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask,
1276 std::list<std::pair<int, int> > &lines2update,
1277 int max_box_nchan)
1278 throw (casa::AipsError)
1279{
1280 try {
1281 if (!lines2update.size()) return; // ignore an empty list
1282
1283 // current line
1284 std::list<std::pair<int,int> >::iterator li=lines2update.begin();
1285 // strongest line
1286 std::list<std::pair<int,int> >::iterator strongli=lines2update.begin();
1287 // the flux (absolute value) of the strongest line
1288 Float peak_flux=-1; // negative value - a flag showing uninitialized
1289 // value
1290 // the algorithm below relies on the list being ordered
1291 Float tmp_flux=-1; // a temporary peak
1292 for (RunningBox running_box(spectrum,temp_mask,edge,max_box_nchan);
1293 running_box.haveMore(); running_box.next()) {
1294
1295 if (li==lines2update.end()) break; // no more lines
1296 const int ch=running_box.getChannel();
1297 if (ch>=li->first && ch<li->second)
1298 if (temp_mask[ch] && tmp_flux<fabs(running_box.aboveMean()))
1299 tmp_flux=fabs(running_box.aboveMean());
1300 if (ch==li->second-1) {
1301 if (peak_flux<tmp_flux) { // if peak_flux=-1, this condition
1302 peak_flux=tmp_flux; // will be satisfied
1303 strongli=li;
1304 }
1305 ++li;
1306 tmp_flux=-1;
1307 }
1308 }
1309 std::list<std::pair<int,int> > res;
1310 res.splice(res.end(),lines2update,strongli);
1311 lines2update.clear();
1312 lines2update.splice(lines2update.end(),res);
1313 }
1314 catch (const AipsError &ae) {
1315 throw;
1316 }
1317 catch (const exception &ex) {
1318 throw AipsError(String("STLineFinder::keepStrongestOnly - STL error: ")+ex.what());
1319 }
1320
1321}
1322
1323//
1324///////////////////////////////////////////////////////////////////////////////
1325
1326
1327///////////////////////////////////////////////////////////////////////////////
1328//
1329// LFLineListOperations - a class incapsulating operations with line lists
1330// The LF prefix stands for Line Finder
1331//
1332
1333// concatenate two lists preserving the order. If two lines appear to
1334// be adjacent, they are joined into the new one
1335void LFLineListOperations::addNewSearchResult(const std::list<pair<int, int> > &newlines,
1336 std::list<std::pair<int, int> > &lines_list)
1337 throw(AipsError)
1338{
1339 try {
1340 for (std::list<pair<int,int> >::const_iterator cli=newlines.begin();
1341 cli!=newlines.end();++cli) {
1342
1343 // the first item, which has a non-void intersection or touches
1344 // the new line
1345 std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(),
1346 lines_list.end(), IntersectsWith(*cli));
1347 // the last such item
1348 std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg,
1349 lines_list.end(), not1(IntersectsWith(*cli)));
1350
1351 // extract all lines which intersect or touch a new one into
1352 // a temporary buffer. This may invalidate the iterators
1353 // line_buffer may be empty, if no lines intersects with a new
1354 // one.
1355 std::list<pair<int,int> > lines_buffer;
1356 lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end);
1357
1358 // build a union of all intersecting lines
1359 pair<int,int> union_line=for_each(lines_buffer.begin(),
1360 lines_buffer.end(),BuildUnion(*cli)).result();
1361
1362 // search for a right place for the new line (union_line) and add
1363 std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(),
1364 lines_list.end(), LaterThan(union_line));
1365 lines_list.insert(pos2insert,union_line);
1366 }
1367 }
1368 catch (const AipsError &ae) {
1369 throw;
1370 }
1371 catch (const exception &ex) {
1372 throw AipsError(String("LFLineListOperations::addNewSearchResult - STL error: ")+ex.what());
1373 }
1374}
1375
1376// extend all line ranges to the point where a value stored in the
1377// specified vector changes (e.g. value-mean change its sign)
1378// This operation is necessary to include line wings, which are below
1379// the detection threshold. If lines becomes adjacent, they are
1380// merged together. Any masked channel stops the extension
1381void LFLineListOperations::searchForWings(std::list<std::pair<int, int> > &newlines,
1382 const casa::Vector<casa::Int> &signs,
1383 const casa::Vector<casa::Bool> &mask,
1384 const std::pair<int,int> &edge) throw(casa::AipsError)
1385{
1386 try {
1387 for (std::list<pair<int,int> >::iterator li=newlines.begin();
1388 li!=newlines.end();++li) {
1389 // update the left hand side
1390 for (int n=li->first-1;n>=edge.first;--n) {
1391 if (!mask[n]) break;
1392 if (signs[n]==signs[li->first] && signs[li->first])
1393 li->first=n;
1394 else break;
1395 }
1396 // update the right hand side
1397 for (int n=li->second;n<edge.second;++n) {
1398 if (!mask[n]) break;
1399 if (signs[n]==signs[li->second-1] && signs[li->second-1])
1400 li->second=n;
1401 else break;
1402 }
1403 }
1404 // need to search for possible mergers.
1405 std::list<std::pair<int, int> > result_buffer;
1406 addNewSearchResult(newlines,result_buffer);
1407 newlines.clear();
1408 newlines.splice(newlines.end(),result_buffer);
1409 }
1410 catch (const AipsError &ae) {
1411 throw;
1412 }
1413 catch (const exception &ex) {
1414 throw AipsError(String("LFLineListOperations::extendLines - STL error: ")+ex.what());
1415 }
1416}
1417
1418//
1419///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.