source: trunk/src/STLineFinder.cpp@ 3030

Last change on this file since 3030 was 3029, checked in by Kana Sugimoto, 10 years ago

New Development: Yes

JIRA Issue: Yes (CAS-6929)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): asap as a whole

Description: committing Darrell's changes to make asap work with merged casacore.


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