source: branches/casa-prerelease/pre-asap/src/STLineFinder.cpp@ 2061

Last change on this file since 2061 was 2012, checked in by WataruKawasaki, 14 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2373, CAS-2620)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: For Scantable::polyBaseline(), parameters and return value have been changed.

Test Programs:

Put in Release Notes: Yes

Module(s): sdbaseline, sd.linefinder

Description: (1) CAS-2373-related:

(1.1) Modified Scantable::polyBaseline() to have the row-based loop inside.

Now it fits and subtracts baseline for all rows and also output info
about the fitting result to logger and text file, while in the
previous version this method just did baseline fitting/subtraction
for one row only and had to be put inside a row-based loop at the
python side ("poly_baseline()" in scantable.py) and result output had
also to be controlled at the python side. Using a test data from NRO
45m telescope (348,000 rows, 512 channels), the processing time of
scantable.poly_baseline() has reduced from 130 minutes to 5-6 minutes.

(1.2) For accelerating the another polynomial fitting function, namely

scantable.auto_poly_baseline(), added a method
Scantable::autoPolyBaseline(). This basically does the same thing
with Scantable::polyBaseline(), but uses linefinfer also to
automatically flag the line regions.

(1.3) To make linefinder usable in Scantable class, added a method

linefinder.setdata(). This makes it possible to apply linefinder
for a float-list data given as a parameter, without setting scantable,
while it was indispensable to set scantable to use linefinger previously.

(2) CAS-2620-related:

Added Scantable::cubicSplineBaseline() and autoCubicSplineBaseline() for
fit baseline using the cubic spline function. Parameters include npiece
(number of polynomial pieces), clipthresh (clipping threshold), and
clipniter (maximum iteration number).



  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 52.6 KB
Line 
1//#---------------------------------------------------------------------------
2//# STLineFinder.cc: A class for automated spectral line search
3//#--------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id: STLineFinder.cpp 2012 2011-02-25 05:51:50Z WataruKawasaki $
30//#---------------------------------------------------------------------------
31
32
33// ASAP
34#include "STLineFinder.h"
35#include "STFitter.h"
36#include "IndexedCompare.h"
37
38// STL
39#include <functional>
40#include <algorithm>
41#include <iostream>
42#include <fstream>
43
44using namespace asap;
45using namespace casa;
46using namespace std;
47
48namespace asap {
49
50///////////////////////////////////////////////////////////////////////////////
51//
52// RunningBox - a running box calculator. This class implements
53// iterations over the specified spectrum and calculates
54// running box filter statistics.
55//
56
57class RunningBox {
58 // The input data to work with. Use reference symantics to avoid
59 // an unnecessary copying
60 const casa::Vector<casa::Float> &spectrum; // a buffer for the spectrum
61 const casa::Vector<casa::Bool> &mask; // associated mask
62 const std::pair<int,int> &edge; // start and stop+1 channels
63 // to work with
64
65 // statistics for running box filtering
66 casa::Float sumf; // sum of fluxes
67 casa::Float sumf2; // sum of squares of fluxes
68 casa::Float sumch; // sum of channel numbers (for linear fit)
69 casa::Float sumch2; // sum of squares of channel numbers (for linear fit)
70 casa::Float sumfch; // sum of flux*(channel number) (for linear fit)
71
72 int box_chan_cntr; // actual number of channels in the box
73 int max_box_nchan; // maximum allowed number of channels in the box
74 // (calculated from boxsize and actual spectrum size)
75 // cache for derivative statistics
76 mutable casa::Bool need2recalculate; // if true, values of the statistics
77 // below are invalid
78 mutable casa::Float linmean; // a value of the linear fit to the
79 // points in the running box
80 mutable casa::Float linvariance; // the same for variance
81 int cur_channel; // the number of the current channel
82 int start_advance; // number of channel from which the box can
83 // be moved (the middle of the box, if there is no
84 // masking)
85public:
86 // set up the object with the references to actual data
87 // as well as the number of channels in the running box
88 RunningBox(const casa::Vector<casa::Float> &in_spectrum,
89 const casa::Vector<casa::Bool> &in_mask,
90 const std::pair<int,int> &in_edge,
91 int in_max_box_nchan) throw(AipsError);
92
93 // access to the statistics
94 const casa::Float& getLinMean() const throw(AipsError);
95 const casa::Float& getLinVariance() const throw(AipsError);
96 const casa::Float aboveMean() const throw(AipsError);
97 int getChannel() const throw();
98
99 // actual number of channels in the box (max_box_nchan, if no channels
100 // are masked)
101 int getNumberOfBoxPoints() const throw();
102
103 // next channel
104 void next() throw(AipsError);
105
106 // checking whether there are still elements
107 casa::Bool haveMore() const throw();
108
109 // go to start
110 void rewind() throw(AipsError);
111
112protected:
113 // supplementary function to control running mean/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
507const casa::Float RunningBox::aboveMean() const throw(AipsError)
508{
509 DebugAssert(cur_channel<edge.second, AipsError);
510 if (need2recalculate) updateDerivativeStatistics();
511 return spectrum[cur_channel]-linmean;
512}
513
514int RunningBox::getChannel() const throw()
515{
516 return cur_channel;
517}
518
519// actual number of channels in the box (max_box_nchan, if no channels
520// are masked)
521int RunningBox::getNumberOfBoxPoints() const throw()
522{
523 return box_chan_cntr;
524}
525
526// supplementary function to control running mean/median calculations.
527// It adds a specified channel to the running box and
528// removes (ch-max_box_nchan+1)'th channel from there
529// Channels, for which the mask is false or index is beyond the
530// allowed range, are ignored
531void RunningBox::advanceRunningBox(int ch) throw(AipsError)
532{
533 if (ch>=edge.first && ch<edge.second)
534 if (mask[ch]) { // ch is a valid channel
535 ++box_chan_cntr;
536 sumf+=spectrum[ch];
537 sumf2+=square(spectrum[ch]);
538 sumch+=Float(ch);
539 sumch2+=square(Float(ch));
540 sumfch+=spectrum[ch]*Float(ch);
541 need2recalculate=True;
542 }
543 int ch2remove=ch-max_box_nchan;
544 if (ch2remove>=edge.first && ch2remove<edge.second)
545 if (mask[ch2remove]) { // ch2remove is a valid channel
546 --box_chan_cntr;
547 sumf-=spectrum[ch2remove];
548 sumf2-=square(spectrum[ch2remove]);
549 sumch-=Float(ch2remove);
550 sumch2-=square(Float(ch2remove));
551 sumfch-=spectrum[ch2remove]*Float(ch2remove);
552 need2recalculate=True;
553 }
554}
555
556// next channel
557void RunningBox::next() throw(AipsError)
558{
559 AlwaysAssert(cur_channel<edge.second,AipsError);
560 ++cur_channel;
561 if (cur_channel+max_box_nchan/2<edge.second && cur_channel>=start_advance)
562 advanceRunningBox(cur_channel+max_box_nchan/2); // update statistics
563}
564
565// checking whether there are still elements
566casa::Bool RunningBox::haveMore() const throw()
567{
568 return cur_channel<edge.second;
569}
570
571// calculate derivative statistics. This function is const, because
572// it updates the cache only
573void RunningBox::updateDerivativeStatistics() const throw(AipsError)
574{
575 AlwaysAssert(box_chan_cntr, AipsError);
576
577 Float mean=sumf/Float(box_chan_cntr);
578
579 // linear LSF formulae
580 Float meanch=sumch/Float(box_chan_cntr);
581 Float meanch2=sumch2/Float(box_chan_cntr);
582 if (meanch==meanch2 || box_chan_cntr<3) {
583 // vertical line in the spectrum, can't calculate linmean and linvariance
584 linmean=0.;
585 linvariance=0.;
586 } else {
587 Float coeff=(sumfch/Float(box_chan_cntr)-meanch*mean)/
588 (meanch2-square(meanch));
589 linmean=coeff*(Float(cur_channel)-meanch)+mean;
590 linvariance=sumf2/Float(box_chan_cntr)-square(mean)-
591 square(coeff)*(meanch2-square(meanch));
592 if (linvariance<0.) {
593 // this shouldn't happen normally, but could be due to round-off error
594 linvariance = 0;
595 } else {
596 linvariance = sqrt(linvariance);
597 }
598 }
599 need2recalculate=False;
600}
601
602
603//
604///////////////////////////////////////////////////////////////////////////////
605
606///////////////////////////////////////////////////////////////////////////////
607//
608// LFAboveThreshold - a running mean/median algorithm for line detection
609//
610//
611
612
613// set up the detection criterion
614LFAboveThreshold::LFAboveThreshold(std::list<pair<int,int> > &in_lines,
615 int in_min_nchan,
616 casa::Float in_threshold,
617 bool use_median,
618 int noise_sample_size) throw() :
619 min_nchan(in_min_nchan), threshold(in_threshold),
620 lines(in_lines), running_box(NULL), itsUseMedian(use_median),
621 itsNoiseSampleSize(noise_sample_size) {}
622
623LFAboveThreshold::~LFAboveThreshold() throw()
624{
625 if (running_box!=NULL) delete running_box;
626}
627
628// replace the detection criterion
629void LFAboveThreshold::setCriterion(int in_min_nchan, casa::Float in_threshold)
630 throw()
631{
632 min_nchan=in_min_nchan;
633 threshold=in_threshold;
634}
635
636// get the sign of runningBox->aboveMean(). The RunningBox pointer
637// should be defined
638casa::Int LFAboveThreshold::getAboveMeanSign() const throw()
639{
640 const Float buf=running_box->aboveMean();
641 if (buf>0) return 1;
642 if (buf<0) return -1;
643 return 0;
644}
645
646
647// process a channel: update cur_line and is_detected before and
648// add a new line to the list, if necessary
649void LFAboveThreshold::processChannel(Bool detect,
650 const casa::Vector<casa::Bool> &mask) throw(casa::AipsError)
651{
652 try {
653 if (is_detected_before) {
654 // we have to check that the current detection has the
655 // same sign of running_box->aboveMean
656 // otherwise it could be a spurious detection
657 if (last_sign && last_sign!=getAboveMeanSign())
658 detect=False;
659 }
660 if (detect) {
661 last_sign=getAboveMeanSign();
662 if (is_detected_before)
663 cur_line.second=running_box->getChannel()+1;
664 else {
665 is_detected_before=True;
666 cur_line.first=running_box->getChannel();
667 cur_line.second=running_box->getChannel()+1;
668 }
669 } else processCurLine(mask);
670 }
671 catch (const AipsError &ae) {
672 throw;
673 }
674 catch (const exception &ex) {
675 throw AipsError(String("LFAboveThreshold::processChannel - STL error: ")+ex.what());
676 }
677}
678
679// process the interval of channels stored in cur_line
680// if it satisfies the criterion, add this interval as a new line
681void LFAboveThreshold::processCurLine(const casa::Vector<casa::Bool> &mask)
682 throw(casa::AipsError)
683{
684 try {
685 if (is_detected_before) {
686 if (cur_line.second-cur_line.first>=min_nchan) {
687 // it was a detection. We need to change the list
688 Bool add_new_line=False;
689 if (lines.size()) {
690 for (int i=lines.back().second;i<cur_line.first;++i)
691 if (mask[i]) { // one valid channel in between
692 // means that we deal with a separate line
693 add_new_line=True;
694 break;
695 }
696 } else add_new_line=True;
697 if (add_new_line)
698 lines.push_back(cur_line);
699 else lines.back().second=cur_line.second;
700 }
701 is_detected_before=False;
702 }
703 }
704 catch (const AipsError &ae) {
705 throw;
706 }
707 catch (const exception &ex) {
708 throw AipsError(String("LFAboveThreshold::processCurLine - STL error: ")+ex.what());
709 }
710}
711
712// return the array with signs of the value-current mean
713// An element is +1 if value>mean, -1 if less, 0 if equal.
714// This array is updated each time the findLines method is called and
715// is used to search the line wings
716const casa::Vector<Int>& LFAboveThreshold::getSigns() const throw()
717{
718 return signs;
719}
720
721// find spectral lines and add them into list
722void LFAboveThreshold::findLines(const casa::Vector<casa::Float> &spectrum,
723 const casa::Vector<casa::Bool> &mask,
724 const std::pair<int,int> &edge,
725 int max_box_nchan)
726 throw(casa::AipsError)
727{
728 const int minboxnchan=4;
729 try {
730
731 if (running_box!=NULL) delete running_box;
732 running_box=new RunningBox(spectrum,mask,edge,max_box_nchan);
733
734 // determine the off-line variance first
735 // an assumption made: lines occupy a small part of the spectrum
736
737 const size_t noiseSampleSize = itsNoiseSampleSize<0 ? size_t(edge.second-edge.first) :
738 std::min(size_t(itsNoiseSampleSize), size_t(edge.second-edge.first));
739 DebugAssert(noiseSampleSize,AipsError);
740 const bool globalNoise = (size_t(edge.second - edge.first) == noiseSampleSize);
741 LFNoiseEstimator ne(noiseSampleSize);
742
743 for (;running_box->haveMore();running_box->next()) {
744 ne.add(running_box->getLinVariance());
745 if (ne.filledToCapacity()) {
746 break;
747 }
748 }
749
750 Float offline_variance = -1; // just a flag that it is unset
751
752 if (globalNoise) {
753 offline_variance = itsUseMedian ? ne.median() : ne.meanLowest80Percent();
754 }
755
756 // actual search algorithm
757 is_detected_before=False;
758
759 // initiate the signs array
760 signs.resize(spectrum.nelements());
761 signs=Vector<Int>(spectrum.nelements(),0);
762
763 //ofstream os("dbg.dat");
764 for (running_box->rewind();running_box->haveMore();
765 running_box->next()) {
766 const int ch=running_box->getChannel();
767 if (!globalNoise) {
768 // add a next point for a local noise estimate
769 ne.add(running_box->getLinVariance());
770 }
771 if (running_box->getNumberOfBoxPoints()>=minboxnchan) {
772 if (!globalNoise) {
773 offline_variance = itsUseMedian ? ne.median() : ne.meanLowest80Percent();
774 }
775 AlwaysAssert(offline_variance>0.,AipsError);
776 processChannel(mask[ch] && (fabs(running_box->aboveMean()) >=
777 threshold*offline_variance), mask);
778 } else processCurLine(mask); // just finish what was accumulated before
779
780 signs[ch]=getAboveMeanSign();
781 //os<<ch<<" "<<spectrum[ch]<<" "<<fabs(running_box->aboveMean())<<" "<<
782 //threshold*offline_variance<<endl;
783 }
784 if (lines.size())
785 searchForWings(lines,signs,mask,edge);
786 }
787 catch (const AipsError &ae) {
788 throw;
789 }
790 catch (const exception &ex) {
791 throw AipsError(String("LFAboveThreshold::findLines - STL error: ")+ex.what());
792 }
793}
794
795//
796///////////////////////////////////////////////////////////////////////////////
797
798///////////////////////////////////////////////////////////////////////////////
799//
800// LFLineListOperations::IntersectsWith - An auxiliary object function
801// to test whether two lines have a non-void intersection
802//
803
804
805// line1 - range of the first line: start channel and stop+1
806LFLineListOperations::IntersectsWith::IntersectsWith(const std::pair<int,int> &in_line1) :
807 line1(in_line1) {}
808
809
810// return true if line2 intersects with line1 with at least one
811// common channel, and false otherwise
812// line2 - range of the second line: start channel and stop+1
813bool LFLineListOperations::IntersectsWith::operator()(const std::pair<int,int> &line2)
814 const throw()
815{
816 if (line2.second<line1.first) return false; // line2 is at lower channels
817 if (line2.first>line1.second) return false; // line2 is at upper channels
818 return true; // line2 has an intersection or is adjacent to line1
819}
820
821//
822///////////////////////////////////////////////////////////////////////////////
823
824///////////////////////////////////////////////////////////////////////////////
825//
826// LFLineListOperations::BuildUnion - An auxiliary object function to build a union
827// of several lines to account for a possibility of merging the nearby lines
828//
829
830// set an initial line (can be a first line in the sequence)
831LFLineListOperations::BuildUnion::BuildUnion(const std::pair<int,int> &line1) :
832 temp_line(line1) {}
833
834// update temp_line with a union of temp_line and new_line
835// provided there is no gap between the lines
836void LFLineListOperations::BuildUnion::operator()(const std::pair<int,int> &new_line)
837 throw()
838{
839 if (new_line.first<temp_line.first) temp_line.first=new_line.first;
840 if (new_line.second>temp_line.second) temp_line.second=new_line.second;
841}
842
843// return the result (temp_line)
844const std::pair<int,int>& LFLineListOperations::BuildUnion::result() const throw()
845{
846 return temp_line;
847}
848
849//
850///////////////////////////////////////////////////////////////////////////////
851
852///////////////////////////////////////////////////////////////////////////////
853//
854// LFLineListOperations::LaterThan - An auxiliary object function to test whether a
855// specified line is at lower spectral channels (to preserve the order in
856// the line list)
857//
858
859// setup the line to compare with
860LFLineListOperations::LaterThan::LaterThan(const std::pair<int,int> &in_line1) :
861 line1(in_line1) {}
862
863// return true if line2 should be placed later than line1
864// in the ordered list (so, it is at greater channel numbers)
865bool LFLineListOperations::LaterThan::operator()(const std::pair<int,int> &line2)
866 const throw()
867{
868 if (line2.second<line1.first) return false; // line2 is at lower channels
869 if (line2.first>line1.second) return true; // line2 is at upper channels
870
871 // line2 intersects with line1. We should have no such situation in
872 // practice
873 return line2.first>line1.first;
874}
875
876//
877///////////////////////////////////////////////////////////////////////////////
878
879
880///////////////////////////////////////////////////////////////////////////////
881//
882// STLineFinder - a class for automated spectral line search
883//
884//
885
886STLineFinder::STLineFinder() throw() : edge(0,0)
887{
888 setOptions();
889}
890
891// set the parameters controlling algorithm
892// in_threshold a single channel threshold default is sqrt(3), which
893// means together with 3 minimum channels at least 3 sigma
894// detection criterion
895// For bad baseline shape, in_threshold may need to be
896// increased
897// in_min_nchan minimum number of channels above the threshold to report
898// a detection, default is 3
899// in_avg_limit perform the averaging of no more than in_avg_limit
900// adjacent channels to search for broad lines
901// Default is 8, but for a bad baseline shape this
902// parameter should be decreased (may be even down to a
903// minimum of 1 to disable this option) to avoid
904// confusing of baseline undulations with a real line.
905// Setting a very large value doesn't usually provide
906// valid detections.
907// in_box_size the box size for running mean/median calculation. Default is
908// 1./5. of the whole spectrum size
909// in_noise_box the box size for off-line noise estimation (if working with
910// local noise. Negative value means use global noise estimate
911// Default is -1 (i.e. estimate using the whole spectrum)
912// in_median true if median statistics is used as opposed to average of
913// the lowest 80% of deviations (default)
914void STLineFinder::setOptions(const casa::Float &in_threshold,
915 const casa::Int &in_min_nchan,
916 const casa::Int &in_avg_limit,
917 const casa::Float &in_box_size,
918 const casa::Float &in_noise_box,
919 const casa::Bool &in_median) throw()
920{
921 threshold=in_threshold;
922 min_nchan=in_min_nchan;
923 avg_limit=in_avg_limit;
924 box_size=in_box_size;
925 itsNoiseBox = in_noise_box;
926 itsUseMedian = in_median;
927
928 useScantable = true;
929}
930
931STLineFinder::~STLineFinder() throw(AipsError) {}
932
933// set scan to work with (in_scan parameter)
934void STLineFinder::setScan(const ScantableWrapper &in_scan) throw(AipsError)
935{
936 scan=in_scan.getCP();
937 AlwaysAssert(!scan.null(),AipsError);
938}
939
940// set spectrum data to work with. this is a method to allow linefinder work
941// without setting scantable for the purpose of using linefinder inside some
942// method in scantable class. (Dec 22, 2010 by W.Kawasaki)
943void STLineFinder::setData(const std::vector<float> &in_spectrum)
944{
945 spectrum = Vector<Float>(in_spectrum);
946 useScantable = false;
947}
948
949// search for spectral lines. Number of lines found is returned
950// in_edge and in_mask control channel rejection for a given row
951// if in_edge has zero length, all channels chosen by mask will be used
952// if in_edge has one element only, it represents the number of
953// channels to drop from both sides of the spectrum
954// in_edge is introduced for convinience, although all functionality
955// can be achieved using a spectrum mask only
956int STLineFinder::findLines(const std::vector<bool> &in_mask,
957 const std::vector<int> &in_edge,
958 const casa::uInt &whichRow) throw(casa::AipsError)
959{
960 if (useScantable && scan.null())
961 throw AipsError("STLineFinder::findLines - a scan should be set first,"
962 " use set_scan");
963
964 uInt nchan = useScantable ? scan->nchan(scan->getIF(whichRow)) : spectrum.nelements();
965 // set up mask and edge rejection
966 // no mask given...
967 if (in_mask.size() == 0) {
968 mask = Vector<Bool>(nchan,True);
969 } else {
970 // use provided mask
971 mask=Vector<Bool>(in_mask);
972 }
973 if (mask.nelements()!=nchan)
974 throw AipsError("STLineFinder::findLines - in_scan and in_mask, or in_spectrum "
975 "and in_mask have different number of spectral channels.");
976
977 // taking flagged channels into account
978 if (useScantable) {
979 vector<bool> flaggedChannels = scan->getMask(whichRow);
980 if (flaggedChannels.size()) {
981 // there is a mask set for this row
982 if (flaggedChannels.size() != mask.nelements()) {
983 throw AipsError("STLineFinder::findLines - internal inconsistency: number of "
984 "mask elements do not match the number of channels");
985 }
986 for (size_t ch = 0; ch<mask.nelements(); ++ch) {
987 mask[ch] &= flaggedChannels[ch];
988 }
989 }
990 }
991
992 // number of elements in in_edge
993 if (in_edge.size()>2)
994 throw AipsError("STLineFinder::findLines - the length of the in_edge parameter"
995 "should not exceed 2");
996 if (!in_edge.size()) {
997 // all spectra, no rejection
998 edge.first=0;
999 edge.second=nchan;
1000 } else {
1001 edge.first=in_edge[0];
1002 if (edge.first<0)
1003 throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
1004 "number of channels to drop");
1005 if (edge.first>=int(nchan))
1006 throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
1007 if (in_edge.size()==2) {
1008 edge.second=in_edge[1];
1009 if (edge.second<0)
1010 throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
1011 "number of channels to drop");
1012 edge.second=nchan-edge.second;
1013 } else edge.second=nchan-edge.first;
1014 if (edge.second<0 || (edge.first>=edge.second))
1015 throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
1016 }
1017
1018 //
1019 int max_box_nchan=int(nchan*box_size); // number of channels in running
1020 // box
1021 if (max_box_nchan<2)
1022 throw AipsError("STLineFinder::findLines - box_size is too small");
1023
1024 // number of elements in the sample for noise estimate
1025 const int noise_box = itsNoiseBox<0 ? -1 : int(nchan * itsNoiseBox);
1026
1027 if ((noise_box!= -1) and (noise_box<2))
1028 throw AipsError("STLineFinder::findLines - noise_box is supposed to be at least 2 elements");
1029
1030 if (useScantable) {
1031 spectrum.resize();
1032 spectrum = Vector<Float>(scan->getSpectrum(whichRow));
1033 }
1034
1035 lines.resize(0); // search from the scratch
1036 last_row_used=whichRow;
1037 Vector<Bool> temp_mask(mask);
1038
1039 Bool first_pass=True;
1040 Int avg_factor=1; // this number of adjacent channels is averaged together
1041 // the total number of the channels is not altered
1042 // instead, min_nchan is also scaled
1043 // it helps to search for broad lines
1044 Vector<Int> signs; // a buffer for signs of the value - mean quantity
1045 // see LFAboveThreshold for details
1046 // We need only signs resulted from last iteration
1047 // because all previous values may be corrupted by the
1048 // presence of spectral lines
1049 while (true) {
1050 // a buffer for new lines found at this iteration
1051 std::list<pair<int,int> > new_lines;
1052
1053 try {
1054 // line find algorithm
1055 LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold, itsUseMedian,noise_box);
1056 lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan);
1057 signs.resize(lfalg.getSigns().nelements());
1058 signs=lfalg.getSigns();
1059 first_pass=False;
1060 if (!new_lines.size())
1061 throw AipsError("spurious"); // nothing new - use the same
1062 // code as for a real exception
1063 }
1064 catch(const AipsError &ae) {
1065 if (first_pass) throw;
1066 // nothing new - proceed to the next step of averaging, if any
1067 // (to search for broad lines)
1068 if (avg_factor>=avg_limit) break; // averaging up to avg_limit
1069 // adjacent channels,
1070 // stop after that
1071 avg_factor*=2; // twice as more averaging
1072 subtractBaseline(temp_mask,9);
1073 averageAdjacentChannels(temp_mask,avg_factor);
1074 continue;
1075 }
1076 keepStrongestOnly(temp_mask,new_lines,max_box_nchan);
1077 // update the list (lines) merging intervals, if necessary
1078 addNewSearchResult(new_lines,lines);
1079 // get a new mask
1080 temp_mask=getMask();
1081 }
1082
1083 // an additional search for wings because in the presence of very strong
1084 // lines temporary mean used at each iteration will be higher than
1085 // the true mean
1086
1087 if (lines.size())
1088 LFLineListOperations::searchForWings(lines,signs,mask,edge);
1089
1090 return int(lines.size());
1091}
1092
1093// auxiliary function to fit and subtract a polynomial from the current
1094// spectrum. It uses the Fitter class. This action is required before
1095// reducing the spectral resolution if the baseline shape is bad
1096void STLineFinder::subtractBaseline(const casa::Vector<casa::Bool> &temp_mask,
1097 const casa::Int &order) throw(casa::AipsError)
1098{
1099 AlwaysAssert(spectrum.nelements(),AipsError);
1100 // use the fact that temp_mask excludes channels rejected at the edge
1101 Fitter sdf;
1102 std::vector<float> absc(spectrum.nelements());
1103 for (unsigned int i=0;i<absc.size();++i)
1104 absc[i]=float(i)/float(spectrum.nelements());
1105 std::vector<float> spec;
1106 spectrum.tovector(spec);
1107 std::vector<bool> std_mask;
1108 temp_mask.tovector(std_mask);
1109 sdf.setData(absc,spec,std_mask);
1110 sdf.setExpression("poly",order);
1111 if (!sdf.fit()) return; // fit failed, use old spectrum
1112 spectrum=casa::Vector<casa::Float>(sdf.getResidual());
1113}
1114
1115// auxiliary function to average adjacent channels and update the mask
1116// if at least one channel involved in summation is masked, all
1117// output channels will be masked. This function works with the
1118// spectrum and edge fields of this class, but updates the mask
1119// array specified, rather than the field of this class
1120// boxsize - a number of adjacent channels to average
1121void STLineFinder::averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update,
1122 const casa::Int &boxsize)
1123 throw(casa::AipsError)
1124{
1125 DebugAssert(mask2update.nelements()==spectrum.nelements(), AipsError);
1126 DebugAssert(boxsize!=0,AipsError);
1127
1128 for (int n=edge.first;n<edge.second;n+=boxsize) {
1129 DebugAssert(n<spectrum.nelements(),AipsError);
1130 int nboxch=0; // number of channels currently in the box
1131 Float mean=0; // buffer for mean calculations
1132 for (int k=n;k<n+boxsize && k<edge.second;++k)
1133 if (mask2update[k]) { // k is a valid channel
1134 mean+=spectrum[k];
1135 ++nboxch;
1136 }
1137 if (nboxch<boxsize) // mask these channels
1138 for (int k=n;k<n+boxsize && k<edge.second;++k)
1139 mask2update[k]=False;
1140 else {
1141 mean/=Float(boxsize);
1142 for (int k=n;k<n+boxsize && k<edge.second;++k)
1143 spectrum[k]=mean;
1144 }
1145 }
1146}
1147
1148
1149// get the mask to mask out all lines that have been found (default)
1150// if invert=true, only channels belong to lines will be unmasked
1151// Note: all channels originally masked by the input mask (in_mask
1152// in setScan) or dropped out by the edge parameter (in_edge
1153// in setScan) are still excluded regardless on the invert option
1154std::vector<bool> STLineFinder::getMask(bool invert)
1155 const throw(casa::AipsError)
1156{
1157 try {
1158 if (useScantable) {
1159 if (scan.null())
1160 throw AipsError("STLineFinder::getMask - a scan should be set first,"
1161 " use set_scan followed by find_lines");
1162 DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
1163 }
1164 /*
1165 if (!lines.size())
1166 throw AipsError("STLineFinder::getMask - one have to search for "
1167 "lines first, use find_lines");
1168 */
1169 std::vector<bool> res_mask(mask.nelements());
1170 // iterator through lines
1171 std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
1172 for (int ch=0;ch<int(res_mask.size());++ch) {
1173 if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
1174 else if (!mask[ch]) res_mask[ch]=false;
1175 else {
1176 res_mask[ch]=!invert; // no line by default
1177 if (cli!=lines.end())
1178 if (ch>=cli->first && ch<cli->second)
1179 res_mask[ch]=invert; // this is a line
1180 }
1181 if (cli!=lines.end())
1182 if (ch>=cli->second)
1183 ++cli; // next line in the list
1184 }
1185 return res_mask;
1186 }
1187 catch (const AipsError &ae) {
1188 throw;
1189 }
1190 catch (const exception &ex) {
1191 throw AipsError(String("STLineFinder::getMask - STL error: ")+ex.what());
1192 }
1193}
1194
1195// get range for all lines found. The same units as used in the scan
1196// will be returned (e.g. velocity instead of channels).
1197std::vector<double> STLineFinder::getLineRanges()
1198 const throw(casa::AipsError)
1199{
1200 std::vector<double> vel;
1201 if (useScantable) {
1202 // convert to required abscissa units
1203 vel = scan->getAbcissa(last_row_used);
1204 } else {
1205 for (int i = 0; i < spectrum.nelements(); ++i)
1206 vel.push_back((double)i);
1207 }
1208 std::vector<int> ranges=getLineRangesInChannels();
1209 std::vector<double> res(ranges.size());
1210
1211 std::vector<int>::const_iterator cri=ranges.begin();
1212 std::vector<double>::iterator outi=res.begin();
1213 for (;cri!=ranges.end() && outi!=res.end();++cri,++outi)
1214 if (uInt(*cri)>=vel.size())
1215 throw AipsError("STLineFinder::getLineRanges - getAbcissa provided less channels than reqired");
1216 else *outi=vel[*cri];
1217 return res;
1218}
1219
1220// The same as getLineRanges, but channels are always used to specify
1221// the range
1222std::vector<int> STLineFinder::getLineRangesInChannels()
1223 const throw(casa::AipsError)
1224{
1225 try {
1226 if (useScantable) {
1227 if (scan.null())
1228 throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first,"
1229 " use set_scan followed by find_lines");
1230 DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
1231 }
1232
1233 if (!lines.size())
1234 throw AipsError("STLineFinder::getLineRangesInChannels - one have to search for "
1235 "lines first, use find_lines");
1236
1237 std::vector<int> res(2*lines.size());
1238 // iterator through lines & result
1239 std::list<std::pair<int,int> >::const_iterator cli = lines.begin();
1240 std::vector<int>::iterator ri = res.begin();
1241 for (; cli != lines.end() && ri != res.end(); ++cli,++ri) {
1242 *ri = cli->first;
1243 if (++ri != res.end())
1244 *ri = cli->second - 1;
1245 }
1246 return res;
1247 } catch (const AipsError &ae) {
1248 throw;
1249 } catch (const exception &ex) {
1250 throw AipsError(String("STLineFinder::getLineRanges - STL error: ") + ex.what());
1251 }
1252}
1253
1254
1255
1256// an auxiliary function to remove all lines from the list, except the
1257// strongest one (by absolute value). If the lines removed are real,
1258// they will be find again at the next iteration. This approach
1259// increases the number of iterations required, but is able to remove
1260// spurious detections likely to occur near strong lines.
1261// Later a better criterion may be implemented, e.g.
1262// taking into consideration the brightness of different lines. Now
1263// use the simplest solution
1264// temp_mask - mask to work with (may be different from original mask as
1265// the lines previously found may be masked)
1266// lines2update - a list of lines to work with
1267// nothing will be done if it is empty
1268// max_box_nchan - channels in the running box for baseline filtering
1269void STLineFinder::keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask,
1270 std::list<std::pair<int, int> > &lines2update,
1271 int max_box_nchan)
1272 throw (casa::AipsError)
1273{
1274 try {
1275 if (!lines2update.size()) return; // ignore an empty list
1276
1277 // current line
1278 std::list<std::pair<int,int> >::iterator li=lines2update.begin();
1279 // strongest line
1280 std::list<std::pair<int,int> >::iterator strongli=lines2update.begin();
1281 // the flux (absolute value) of the strongest line
1282 Float peak_flux=-1; // negative value - a flag showing uninitialized
1283 // value
1284 // the algorithm below relies on the list being ordered
1285 Float tmp_flux=-1; // a temporary peak
1286 for (RunningBox running_box(spectrum,temp_mask,edge,max_box_nchan);
1287 running_box.haveMore(); running_box.next()) {
1288
1289 if (li==lines2update.end()) break; // no more lines
1290 const int ch=running_box.getChannel();
1291 if (ch>=li->first && ch<li->second)
1292 if (temp_mask[ch] && tmp_flux<fabs(running_box.aboveMean()))
1293 tmp_flux=fabs(running_box.aboveMean());
1294 if (ch==li->second-1) {
1295 if (peak_flux<tmp_flux) { // if peak_flux=-1, this condition
1296 peak_flux=tmp_flux; // will be satisfied
1297 strongli=li;
1298 }
1299 ++li;
1300 tmp_flux=-1;
1301 }
1302 }
1303 std::list<std::pair<int,int> > res;
1304 res.splice(res.end(),lines2update,strongli);
1305 lines2update.clear();
1306 lines2update.splice(lines2update.end(),res);
1307 }
1308 catch (const AipsError &ae) {
1309 throw;
1310 }
1311 catch (const exception &ex) {
1312 throw AipsError(String("STLineFinder::keepStrongestOnly - STL error: ")+ex.what());
1313 }
1314
1315}
1316
1317//
1318///////////////////////////////////////////////////////////////////////////////
1319
1320
1321///////////////////////////////////////////////////////////////////////////////
1322//
1323// LFLineListOperations - a class incapsulating operations with line lists
1324// The LF prefix stands for Line Finder
1325//
1326
1327// concatenate two lists preserving the order. If two lines appear to
1328// be adjacent, they are joined into the new one
1329void LFLineListOperations::addNewSearchResult(const std::list<pair<int, int> > &newlines,
1330 std::list<std::pair<int, int> > &lines_list)
1331 throw(AipsError)
1332{
1333 try {
1334 for (std::list<pair<int,int> >::const_iterator cli=newlines.begin();
1335 cli!=newlines.end();++cli) {
1336
1337 // the first item, which has a non-void intersection or touches
1338 // the new line
1339 std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(),
1340 lines_list.end(), IntersectsWith(*cli));
1341 // the last such item
1342 std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg,
1343 lines_list.end(), not1(IntersectsWith(*cli)));
1344
1345 // extract all lines which intersect or touch a new one into
1346 // a temporary buffer. This may invalidate the iterators
1347 // line_buffer may be empty, if no lines intersects with a new
1348 // one.
1349 std::list<pair<int,int> > lines_buffer;
1350 lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end);
1351
1352 // build a union of all intersecting lines
1353 pair<int,int> union_line=for_each(lines_buffer.begin(),
1354 lines_buffer.end(),BuildUnion(*cli)).result();
1355
1356 // search for a right place for the new line (union_line) and add
1357 std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(),
1358 lines_list.end(), LaterThan(union_line));
1359 lines_list.insert(pos2insert,union_line);
1360 }
1361 }
1362 catch (const AipsError &ae) {
1363 throw;
1364 }
1365 catch (const exception &ex) {
1366 throw AipsError(String("LFLineListOperations::addNewSearchResult - STL error: ")+ex.what());
1367 }
1368}
1369
1370// extend all line ranges to the point where a value stored in the
1371// specified vector changes (e.g. value-mean change its sign)
1372// This operation is necessary to include line wings, which are below
1373// the detection threshold. If lines becomes adjacent, they are
1374// merged together. Any masked channel stops the extension
1375void LFLineListOperations::searchForWings(std::list<std::pair<int, int> > &newlines,
1376 const casa::Vector<casa::Int> &signs,
1377 const casa::Vector<casa::Bool> &mask,
1378 const std::pair<int,int> &edge) throw(casa::AipsError)
1379{
1380 try {
1381 for (std::list<pair<int,int> >::iterator li=newlines.begin();
1382 li!=newlines.end();++li) {
1383 // update the left hand side
1384 for (int n=li->first-1;n>=edge.first;--n) {
1385 if (!mask[n]) break;
1386 if (signs[n]==signs[li->first] && signs[li->first])
1387 li->first=n;
1388 else break;
1389 }
1390 // update the right hand side
1391 for (int n=li->second;n<edge.second;++n) {
1392 if (!mask[n]) break;
1393 if (signs[n]==signs[li->second-1] && signs[li->second-1])
1394 li->second=n;
1395 else break;
1396 }
1397 }
1398 // need to search for possible mergers.
1399 std::list<std::pair<int, int> > result_buffer;
1400 addNewSearchResult(newlines,result_buffer);
1401 newlines.clear();
1402 newlines.splice(newlines.end(),result_buffer);
1403 }
1404 catch (const AipsError &ae) {
1405 throw;
1406 }
1407 catch (const exception &ex) {
1408 throw AipsError(String("LFLineListOperations::extendLines - STL error: ")+ex.what());
1409 }
1410}
1411
1412//
1413///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.