source: branches/alma/src/STLineFinder.cpp@ 2521

Last change on this file since 2521 was 1757, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 52.3 KB
RevLine 
[297]1//#---------------------------------------------------------------------------
[881]2//# STLineFinder.cc: A class for automated spectral line search
[297]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//#
[890]29//# $Id: STLineFinder.cpp 1757 2010-06-09 09:03:06Z KanaSugimoto $
[297]30//#---------------------------------------------------------------------------
31
32
33// ASAP
[894]34#include "STLineFinder.h"
35#include "STFitter.h"
[1757]36#include "IndexedCompare.h"
[297]37
38// STL
[343]39#include <functional>
40#include <algorithm>
[297]41#include <iostream>
[351]42#include <fstream>
[297]43
44using namespace asap;
45using namespace casa;
46using namespace std;
47
[344]48namespace asap {
49
[343]50///////////////////////////////////////////////////////////////////////////////
51//
[881]52// RunningBox - a running box calculator. This class implements
[1315]53// iterations over the specified spectrum and calculates
[351]54// running box filter statistics.
[343]55//
56
[351]57class RunningBox {
[331]58 // The input data to work with. Use reference symantics to avoid
[881]59 // an unnecessary copying
[331]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
[881]64
[351]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)
[881]71
[331]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)
[351]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
[996]84 // masking)
[351]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,
[996]90 const std::pair<int,int> &in_edge,
91 int in_max_box_nchan) throw(AipsError);
[881]92
[351]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();
[881]98
[351]99 // actual number of channels in the box (max_box_nchan, if no channels
100 // are masked)
101 int getNumberOfBoxPoints() const throw();
[297]102
[351]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);
[881]111
[351]112protected:
[1757]113 // supplementary function to control running mean/median calculations.
114 // It adds a specified channel to the running box and
[351]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//
[352]135class LFAboveThreshold : protected LFLineListOperations {
[331]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
[996]142 // the detection criterion, to be
143 // a detection
[881]144 casa::Float threshold; // detection threshold - the
[331]145 // minimal signal to noise ratio
[351]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
[551]149 casa::Vector<Int> signs; // An array to store the signs of
150 // the value - current mean
[996]151 // (used to search wings)
[907]152 casa::Int last_sign; // a sign (+1, -1 or 0) of the
153 // last point of the detected line
154 //
[1757]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)
[331]161public:
[351]162
163 // set up the detection criterion
164 LFAboveThreshold(std::list<pair<int,int> > &in_lines,
165 int in_min_nchan = 3,
[1757]166 casa::Float in_threshold = 5,
167 bool use_median = false,
168 int noise_sample_size = -1) throw();
[351]169 virtual ~LFAboveThreshold() throw();
[881]170
[331]171 // replace the detection criterion
172 void setCriterion(int in_min_nchan, casa::Float in_threshold) throw();
[297]173
[551]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
[331]180 // find spectral lines and add them into list
[344]181 // if statholder is not NULL, the accumulate function of it will be
182 // called for each channel to save statistics
[351]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,
[996]186 const casa::Vector<casa::Bool> &mask,
187 const std::pair<int,int> &edge,
188 int max_box_nchan) throw(casa::AipsError);
[351]189
[331]190protected:
[297]191
[331]192 // process a channel: update curline and is_detected before and
193 // add a new line to the list, if necessary using processCurLine()
[351]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);
[297]197
[331]198 // process the interval of channels stored in curline
199 // if it satisfies the criterion, add this interval as a new line
[351]200 void processCurLine(const casa::Vector<casa::Bool> &mask)
201 throw(casa::AipsError);
[924]202
[907]203 // get the sign of runningBox->aboveMean(). The RunningBox pointer
204 // should be defined
205 casa::Int getAboveMeanSign() const throw();
[331]206};
[344]207
208//
209///////////////////////////////////////////////////////////////////////////////
[351]210
[1757]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
[331]278} // namespace asap
[297]279
[344]280///////////////////////////////////////////////////////////////////////////////
[343]281//
[1757]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//
[881]455// RunningBox - a running box calculator. This class implements
[351]456// interations over the specified spectrum and calculates
457// running box filter statistics.
[331]458//
[297]459
[331]460// set up the object with the references to actual data
461// and the number of channels in the running box
[351]462RunningBox::RunningBox(const casa::Vector<casa::Float> &in_spectrum,
463 const casa::Vector<casa::Bool> &in_mask,
[996]464 const std::pair<int,int> &in_edge,
465 int in_max_box_nchan) throw(AipsError) :
[331]466 spectrum(in_spectrum), mask(in_mask), edge(in_edge),
[996]467 max_box_nchan(in_max_box_nchan)
[351]468{
469 rewind();
470}
[331]471
[351]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);
[881]484
485 if (initial_box_ch==edge.second)
[351]486 throw AipsError("RunningBox::rewind - too much channels are masked");
487
488 cur_channel=edge.first;
[881]489 start_advance=initial_box_ch-max_box_nchan/2;
[351]490}
491
492// access to the statistics
493const casa::Float& RunningBox::getLinMean() const throw(AipsError)
[331]494{
[351]495 DebugAssert(cur_channel<edge.second, AipsError);
496 if (need2recalculate) updateDerivativeStatistics();
497 return linmean;
[297]498}
499
[351]500const casa::Float& RunningBox::getLinVariance() const throw(AipsError)
501{
502 DebugAssert(cur_channel<edge.second, AipsError);
503 if (need2recalculate) updateDerivativeStatistics();
504 return linvariance;
505}
[331]506
[351]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
[1757]526// supplementary function to control running mean/median calculations.
527// It adds a specified channel to the running box and
[297]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
[351]531void RunningBox::advanceRunningBox(int ch) throw(AipsError)
[297]532{
533 if (ch>=edge.first && ch<edge.second)
534 if (mask[ch]) { // ch is a valid channel
535 ++box_chan_cntr;
[351]536 sumf+=spectrum[ch];
537 sumf2+=square(spectrum[ch]);
[996]538 sumch+=Float(ch);
539 sumch2+=square(Float(ch));
540 sumfch+=spectrum[ch]*Float(ch);
541 need2recalculate=True;
[297]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;
[351]547 sumf-=spectrum[ch2remove];
[881]548 sumf2-=square(spectrum[ch2remove]);
[996]549 sumch-=Float(ch2remove);
550 sumch2-=square(Float(ch2remove));
551 sumfch-=spectrum[ch2remove]*Float(ch2remove);
552 need2recalculate=True;
[297]553 }
554}
555
[351]556// next channel
557void RunningBox::next() throw(AipsError)
[297]558{
[351]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
[297]563}
564
[351]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);
[881]576
[351]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;
[1757]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 }
[351]598 }
599 need2recalculate=False;
600}
601
602
603//
604///////////////////////////////////////////////////////////////////////////////
605
606///////////////////////////////////////////////////////////////////////////////
607//
[1757]608// LFAboveThreshold - a running mean/median algorithm for line detection
[351]609//
610//
611
612
613// set up the detection criterion
614LFAboveThreshold::LFAboveThreshold(std::list<pair<int,int> > &in_lines,
615 int in_min_nchan,
[1757]616 casa::Float in_threshold,
617 bool use_median,
618 int noise_sample_size) throw() :
[351]619 min_nchan(in_min_nchan), threshold(in_threshold),
[1757]620 lines(in_lines), running_box(NULL), itsUseMedian(use_median),
621 itsNoiseSampleSize(noise_sample_size) {}
[351]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
[907]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}
[351]645
[907]646
[297]647// process a channel: update cur_line and is_detected before and
648// add a new line to the list, if necessary
[351]649void LFAboveThreshold::processChannel(Bool detect,
650 const casa::Vector<casa::Bool> &mask) throw(casa::AipsError)
[297]651{
652 try {
[907]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;
[1315]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);
[297]670 }
671 catch (const AipsError &ae) {
672 throw;
[881]673 }
[297]674 catch (const exception &ex) {
[351]675 throw AipsError(String("LFAboveThreshold::processChannel - STL error: ")+ex.what());
[297]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
[351]681void LFAboveThreshold::processCurLine(const casa::Vector<casa::Bool> &mask)
[331]682 throw(casa::AipsError)
[297]683{
684 try {
[881]685 if (is_detected_before) {
[1315]686 if (cur_line.second-cur_line.first>=min_nchan) {
[996]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);
[881]699 else lines.back().second=cur_line.second;
[996]700 }
701 is_detected_before=False;
[881]702 }
[297]703 }
704 catch (const AipsError &ae) {
705 throw;
[881]706 }
[297]707 catch (const exception &ex) {
[351]708 throw AipsError(String("LFAboveThreshold::processCurLine - STL error: ")+ex.what());
[297]709 }
710}
711
[551]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
[331]721// find spectral lines and add them into list
[351]722void LFAboveThreshold::findLines(const casa::Vector<casa::Float> &spectrum,
[996]723 const casa::Vector<casa::Bool> &mask,
724 const std::pair<int,int> &edge,
725 int max_box_nchan)
[331]726 throw(casa::AipsError)
727{
728 const int minboxnchan=4;
[351]729 try {
[331]730
[351]731 if (running_box!=NULL) delete running_box;
732 running_box=new RunningBox(spectrum,mask,edge,max_box_nchan);
[368]733
734 // determine the off-line variance first
735 // an assumption made: lines occupy a small part of the spectrum
[881]736
[1757]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);
[881]742
[1757]743 for (;running_box->haveMore();running_box->next()) {
744 ne.add(running_box->getLinVariance());
745 if (ne.filledToCapacity()) {
746 break;
747 }
748 }
[881]749
[1757]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 }
[881]755
[351]756 // actual search algorithm
757 is_detected_before=False;
[368]758
[551]759 // initiate the signs array
760 signs.resize(spectrum.nelements());
761 signs=Vector<Int>(spectrum.nelements(),0);
762
[369]763 //ofstream os("dbg.dat");
[368]764 for (running_box->rewind();running_box->haveMore();
765 running_box->next()) {
[351]766 const int ch=running_box->getChannel();
[1757]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);
[996]776 processChannel(mask[ch] && (fabs(running_box->aboveMean()) >=
777 threshold*offline_variance), mask);
[1757]778 } else processCurLine(mask); // just finish what was accumulated before
[907]779
[996]780 signs[ch]=getAboveMeanSign();
[1757]781 //os<<ch<<" "<<spectrum[ch]<<" "<<fabs(running_box->aboveMean())<<" "<<
782 //threshold*offline_variance<<endl;
[351]783 }
[352]784 if (lines.size())
785 searchForWings(lines,signs,mask,edge);
[344]786 }
[351]787 catch (const AipsError &ae) {
788 throw;
[881]789 }
[351]790 catch (const exception &ex) {
791 throw AipsError(String("LFAboveThreshold::findLines - STL error: ")+ex.what());
792 }
[331]793}
794
795//
796///////////////////////////////////////////////////////////////////////////////
797
[343]798///////////////////////////////////////////////////////////////////////////////
799//
[352]800// LFLineListOperations::IntersectsWith - An auxiliary object function
801// to test whether two lines have a non-void intersection
[343]802//
[331]803
[343]804
805// line1 - range of the first line: start channel and stop+1
[352]806LFLineListOperations::IntersectsWith::IntersectsWith(const std::pair<int,int> &in_line1) :
[343]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
[352]813bool LFLineListOperations::IntersectsWith::operator()(const std::pair<int,int> &line2)
[343]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//
[352]826// LFLineListOperations::BuildUnion - An auxiliary object function to build a union
[343]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)
[352]831LFLineListOperations::BuildUnion::BuildUnion(const std::pair<int,int> &line1) :
[343]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
[352]836void LFLineListOperations::BuildUnion::operator()(const std::pair<int,int> &new_line)
[343]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)
[352]844const std::pair<int,int>& LFLineListOperations::BuildUnion::result() const throw()
[343]845{
846 return temp_line;
847}
848
849//
850///////////////////////////////////////////////////////////////////////////////
851
852///////////////////////////////////////////////////////////////////////////////
853//
[352]854// LFLineListOperations::LaterThan - An auxiliary object function to test whether a
[343]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
[352]860LFLineListOperations::LaterThan::LaterThan(const std::pair<int,int> &in_line1) :
[343]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)
[352]865bool LFLineListOperations::LaterThan::operator()(const std::pair<int,int> &line2)
[343]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
[881]870
[343]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//
[881]882// STLineFinder - a class for automated spectral line search
[343]883//
884//
[331]885
[881]886STLineFinder::STLineFinder() throw() : edge(0,0)
[331]887{
[369]888 setOptions();
[331]889}
890
[369]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
[881]901// Default is 8, but for a bad baseline shape this
[369]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.
[881]905// Setting a very large value doesn't usually provide
906// valid detections.
[1757]907// in_box_size the box size for running mean/median calculation. Default is
[369]908// 1./5. of the whole spectrum size
[1757]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)
[881]914void STLineFinder::setOptions(const casa::Float &in_threshold,
[369]915 const casa::Int &in_min_nchan,
[996]916 const casa::Int &in_avg_limit,
[1757]917 const casa::Float &in_box_size,
918 const casa::Float &in_noise_box,
919 const casa::Bool &in_median) throw()
[369]920{
921 threshold=in_threshold;
922 min_nchan=in_min_nchan;
923 avg_limit=in_avg_limit;
924 box_size=in_box_size;
[1757]925 itsNoiseBox = in_noise_box;
926 itsUseMedian = in_median;
[369]927}
928
[881]929STLineFinder::~STLineFinder() throw(AipsError) {}
[331]930
[907]931// set scan to work with (in_scan parameter)
932void STLineFinder::setScan(const ScantableWrapper &in_scan) throw(AipsError)
933{
934 scan=in_scan.getCP();
935 AlwaysAssert(!scan.null(),AipsError);
[924]936
[907]937}
938
939// search for spectral lines. Number of lines found is returned
940// in_edge and in_mask control channel rejection for a given row
[331]941// if in_edge has zero length, all channels chosen by mask will be used
942// if in_edge has one element only, it represents the number of
943// channels to drop from both sides of the spectrum
944// in_edge is introduced for convinience, although all functionality
[881]945// can be achieved using a spectrum mask only
[907]946int STLineFinder::findLines(const std::vector<bool> &in_mask,
[996]947 const std::vector<int> &in_edge,
948 const casa::uInt &whichRow) throw(casa::AipsError)
[331]949{
[907]950 if (scan.null())
951 throw AipsError("STLineFinder::findLines - a scan should be set first,"
952 " use set_scan");
[924]953
954 uInt nchan = scan->nchan(scan->getIF(whichRow));
[907]955 // set up mask and edge rejection
[924]956 // no mask given...
957 if (in_mask.size() == 0) {
958 mask = Vector<Bool>(nchan,True);
959 } else {
960 // use provided mask
961 mask=Vector<Bool>(in_mask);
962 }
963 if (mask.nelements()!=nchan)
[907]964 throw AipsError("STLineFinder::findLines - in_scan and in_mask have different"
965 "number of spectral channels.");
[1757]966
967 // taking flagged channels into account
968 vector<bool> flaggedChannels = scan->getMask(whichRow);
969 if (flaggedChannels.size()) {
970 // there is a mask set for this row
971 if (flaggedChannels.size() != mask.nelements()) {
972 throw AipsError("STLineFinder::findLines - internal inconsistency: number of mask elements do not match the number of channels");
973 }
974 for (size_t ch = 0; ch<mask.nelements(); ++ch) {
975 mask[ch] &= flaggedChannels[ch];
976 }
977 }
978
[907]979 // number of elements in in_edge
980 if (in_edge.size()>2)
981 throw AipsError("STLineFinder::findLines - the length of the in_edge parameter"
[996]982 "should not exceed 2");
[907]983 if (!in_edge.size()) {
[881]984 // all spectra, no rejection
[331]985 edge.first=0;
[996]986 edge.second=nchan;
[907]987 } else {
988 edge.first=in_edge[0];
[996]989 if (edge.first<0)
990 throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
991 "number of channels to drop");
992 if (edge.first>=int(nchan))
993 throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
[907]994 if (in_edge.size()==2) {
[996]995 edge.second=in_edge[1];
996 if (edge.second<0)
997 throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
998 "number of channels to drop");
[924]999 edge.second=nchan-edge.second;
[996]1000 } else edge.second=nchan-edge.first;
[369]1001 if (edge.second<0 || (edge.first>=edge.second))
[996]1002 throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
[881]1003 }
[924]1004
[907]1005 //
[924]1006 int max_box_nchan=int(nchan*box_size); // number of channels in running
[331]1007 // box
1008 if (max_box_nchan<2)
[881]1009 throw AipsError("STLineFinder::findLines - box_size is too small");
[331]1010
[1757]1011 // number of elements in the sample for noise estimate
1012 const int noise_box = itsNoiseBox<0 ? -1 : int(nchan * itsNoiseBox);
1013
1014 if ((noise_box!= -1) and (noise_box<2))
1015 throw AipsError("STLineFinder::findLines - noise_box is supposed to be at least 2 elements");
1016
[881]1017 spectrum.resize();
1018 spectrum = Vector<Float>(scan->getSpectrum(whichRow));
[331]1019
1020 lines.resize(0); // search from the scratch
[370]1021 last_row_used=whichRow;
[331]1022 Vector<Bool> temp_mask(mask);
[351]1023
1024 Bool first_pass=True;
[368]1025 Int avg_factor=1; // this number of adjacent channels is averaged together
1026 // the total number of the channels is not altered
[996]1027 // instead, min_nchan is also scaled
1028 // it helps to search for broad lines
[551]1029 Vector<Int> signs; // a buffer for signs of the value - mean quantity
1030 // see LFAboveThreshold for details
[996]1031 // We need only signs resulted from last iteration
1032 // because all previous values may be corrupted by the
1033 // presence of spectral lines
[344]1034 while (true) {
[351]1035 // a buffer for new lines found at this iteration
[881]1036 std::list<pair<int,int> > new_lines;
[351]1037
1038 try {
[369]1039 // line find algorithm
[1757]1040 LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold, itsUseMedian,noise_box);
[352]1041 lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan);
[996]1042 signs.resize(lfalg.getSigns().nelements());
1043 signs=lfalg.getSigns();
[368]1044 first_pass=False;
1045 if (!new_lines.size())
[996]1046 throw AipsError("spurious"); // nothing new - use the same
1047 // code as for a real exception
[351]1048 }
1049 catch(const AipsError &ae) {
1050 if (first_pass) throw;
[368]1051 // nothing new - proceed to the next step of averaging, if any
[996]1052 // (to search for broad lines)
[1315]1053 if (avg_factor>=avg_limit) break; // averaging up to avg_limit
[996]1054 // adjacent channels,
1055 // stop after that
1056 avg_factor*=2; // twice as more averaging
1057 subtractBaseline(temp_mask,9);
1058 averageAdjacentChannels(temp_mask,avg_factor);
1059 continue;
[1315]1060 }
[368]1061 keepStrongestOnly(temp_mask,new_lines,max_box_nchan);
[343]1062 // update the list (lines) merging intervals, if necessary
[344]1063 addNewSearchResult(new_lines,lines);
1064 // get a new mask
[881]1065 temp_mask=getMask();
[343]1066 }
[881]1067
[551]1068 // an additional search for wings because in the presence of very strong
1069 // lines temporary mean used at each iteration will be higher than
1070 // the true mean
[881]1071
[551]1072 if (lines.size())
1073 LFLineListOperations::searchForWings(lines,signs,mask,edge);
[881]1074
[331]1075 return int(lines.size());
1076}
1077
[369]1078// auxiliary function to fit and subtract a polynomial from the current
[890]1079// spectrum. It uses the Fitter class. This action is required before
[369]1080// reducing the spectral resolution if the baseline shape is bad
[881]1081void STLineFinder::subtractBaseline(const casa::Vector<casa::Bool> &temp_mask,
[369]1082 const casa::Int &order) throw(casa::AipsError)
1083{
1084 AlwaysAssert(spectrum.nelements(),AipsError);
1085 // use the fact that temp_mask excludes channels rejected at the edge
[890]1086 Fitter sdf;
[369]1087 std::vector<float> absc(spectrum.nelements());
[996]1088 for (unsigned int i=0;i<absc.size();++i)
[369]1089 absc[i]=float(i)/float(spectrum.nelements());
1090 std::vector<float> spec;
1091 spectrum.tovector(spec);
1092 std::vector<bool> std_mask;
1093 temp_mask.tovector(std_mask);
1094 sdf.setData(absc,spec,std_mask);
1095 sdf.setExpression("poly",order);
1096 if (!sdf.fit()) return; // fit failed, use old spectrum
[881]1097 spectrum=casa::Vector<casa::Float>(sdf.getResidual());
[369]1098}
1099
[368]1100// auxiliary function to average adjacent channels and update the mask
1101// if at least one channel involved in summation is masked, all
1102// output channels will be masked. This function works with the
1103// spectrum and edge fields of this class, but updates the mask
1104// array specified, rather than the field of this class
1105// boxsize - a number of adjacent channels to average
[881]1106void STLineFinder::averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update,
[368]1107 const casa::Int &boxsize)
1108 throw(casa::AipsError)
1109{
1110 DebugAssert(mask2update.nelements()==spectrum.nelements(), AipsError);
1111 DebugAssert(boxsize!=0,AipsError);
[881]1112
[368]1113 for (int n=edge.first;n<edge.second;n+=boxsize) {
1114 DebugAssert(n<spectrum.nelements(),AipsError);
1115 int nboxch=0; // number of channels currently in the box
1116 Float mean=0; // buffer for mean calculations
1117 for (int k=n;k<n+boxsize && k<edge.second;++k)
1118 if (mask2update[k]) { // k is a valid channel
[996]1119 mean+=spectrum[k];
1120 ++nboxch;
[881]1121 }
[368]1122 if (nboxch<boxsize) // mask these channels
1123 for (int k=n;k<n+boxsize && k<edge.second;++k)
[996]1124 mask2update[k]=False;
[368]1125 else {
1126 mean/=Float(boxsize);
[996]1127 for (int k=n;k<n+boxsize && k<edge.second;++k)
1128 spectrum[k]=mean;
[368]1129 }
1130 }
1131}
[331]1132
[368]1133
[297]1134// get the mask to mask out all lines that have been found (default)
1135// if invert=true, only channels belong to lines will be unmasked
1136// Note: all channels originally masked by the input mask (in_mask
1137// in setScan) or dropped out by the edge parameter (in_edge
1138// in setScan) are still excluded regardless on the invert option
[881]1139std::vector<bool> STLineFinder::getMask(bool invert)
[297]1140 const throw(casa::AipsError)
1141{
1142 try {
1143 if (scan.null())
[881]1144 throw AipsError("STLineFinder::getMask - a scan should be set first,"
[297]1145 " use set_scan followed by find_lines");
[924]1146 DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
[297]1147 /*
1148 if (!lines.size())
[881]1149 throw AipsError("STLineFinder::getMask - one have to search for "
[996]1150 "lines first, use find_lines");
[881]1151 */
[297]1152 std::vector<bool> res_mask(mask.nelements());
1153 // iterator through lines
1154 std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
[1603]1155 for (int ch=0;ch<int(res_mask.size());++ch) {
[297]1156 if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
[996]1157 else if (!mask[ch]) res_mask[ch]=false;
1158 else {
1159 res_mask[ch]=!invert; // no line by default
[1603]1160 if (cli!=lines.end())
1161 if (ch>=cli->first && ch<cli->second)
1162 res_mask[ch]=invert; // this is a line
1163 }
1164 if (cli!=lines.end())
1165 if (ch>=cli->second) {
1166 ++cli; // next line in the list
1167 }
1168 }
[297]1169 return res_mask;
1170 }
1171 catch (const AipsError &ae) {
1172 throw;
[881]1173 }
[297]1174 catch (const exception &ex) {
[881]1175 throw AipsError(String("STLineFinder::getMask - STL error: ")+ex.what());
[297]1176 }
1177}
1178
[370]1179// get range for all lines found. The same units as used in the scan
1180// will be returned (e.g. velocity instead of channels).
[881]1181std::vector<double> STLineFinder::getLineRanges()
[297]1182 const throw(casa::AipsError)
1183{
[370]1184 // convert to required abscissa units
1185 std::vector<double> vel=scan->getAbcissa(last_row_used);
1186 std::vector<int> ranges=getLineRangesInChannels();
1187 std::vector<double> res(ranges.size());
1188
1189 std::vector<int>::const_iterator cri=ranges.begin();
1190 std::vector<double>::iterator outi=res.begin();
1191 for (;cri!=ranges.end() && outi!=res.end();++cri,++outi)
1192 if (uInt(*cri)>=vel.size())
[881]1193 throw AipsError("STLineFinder::getLineRanges - getAbcissa provided less channels than reqired");
[370]1194 else *outi=vel[*cri];
1195 return res;
1196}
1197
1198// The same as getLineRanges, but channels are always used to specify
1199// the range
[881]1200std::vector<int> STLineFinder::getLineRangesInChannels()
[370]1201 const throw(casa::AipsError)
1202{
[297]1203 try {
1204 if (scan.null())
[881]1205 throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first,"
[297]1206 " use set_scan followed by find_lines");
[924]1207 DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
[881]1208
[297]1209 if (!lines.size())
[881]1210 throw AipsError("STLineFinder::getLineRangesInChannels - one have to search for "
[996]1211 "lines first, use find_lines");
[881]1212
[297]1213 std::vector<int> res(2*lines.size());
1214 // iterator through lines & result
1215 std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
1216 std::vector<int>::iterator ri=res.begin();
[881]1217 for (;cli!=lines.end() && ri!=res.end();++cli,++ri) {
[996]1218 *ri=cli->first;
1219 if (++ri!=res.end())
1220 *ri=cli->second-1;
[881]1221 }
[297]1222 return res;
1223 }
1224 catch (const AipsError &ae) {
1225 throw;
[881]1226 }
[297]1227 catch (const exception &ex) {
[881]1228 throw AipsError(String("STLineFinder::getLineRanges - STL error: ")+ex.what());
[297]1229 }
1230}
[331]1231
[370]1232
1233
[368]1234// an auxiliary function to remove all lines from the list, except the
1235// strongest one (by absolute value). If the lines removed are real,
[881]1236// they will be find again at the next iteration. This approach
1237// increases the number of iterations required, but is able to remove
[1315]1238// spurious detections likely to occur near strong lines.
[368]1239// Later a better criterion may be implemented, e.g.
1240// taking into consideration the brightness of different lines. Now
[881]1241// use the simplest solution
[368]1242// temp_mask - mask to work with (may be different from original mask as
1243// the lines previously found may be masked)
1244// lines2update - a list of lines to work with
1245// nothing will be done if it is empty
1246// max_box_nchan - channels in the running box for baseline filtering
[881]1247void STLineFinder::keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask,
[996]1248 std::list<std::pair<int, int> > &lines2update,
1249 int max_box_nchan)
[368]1250 throw (casa::AipsError)
1251{
1252 try {
1253 if (!lines2update.size()) return; // ignore an empty list
1254
1255 // current line
1256 std::list<std::pair<int,int> >::iterator li=lines2update.begin();
1257 // strongest line
1258 std::list<std::pair<int,int> >::iterator strongli=lines2update.begin();
1259 // the flux (absolute value) of the strongest line
1260 Float peak_flux=-1; // negative value - a flag showing uninitialized
1261 // value
1262 // the algorithm below relies on the list being ordered
1263 Float tmp_flux=-1; // a temporary peak
1264 for (RunningBox running_box(spectrum,temp_mask,edge,max_box_nchan);
1265 running_box.haveMore(); running_box.next()) {
1266
1267 if (li==lines2update.end()) break; // no more lines
[996]1268 const int ch=running_box.getChannel();
1269 if (ch>=li->first && ch<li->second)
1270 if (temp_mask[ch] && tmp_flux<fabs(running_box.aboveMean()))
1271 tmp_flux=fabs(running_box.aboveMean());
1272 if (ch==li->second-1) {
1273 if (peak_flux<tmp_flux) { // if peak_flux=-1, this condition
1274 peak_flux=tmp_flux; // will be satisfied
1275 strongli=li;
1276 }
1277 ++li;
1278 tmp_flux=-1;
1279 }
[881]1280 }
[368]1281 std::list<std::pair<int,int> > res;
1282 res.splice(res.end(),lines2update,strongli);
1283 lines2update.clear();
1284 lines2update.splice(lines2update.end(),res);
1285 }
1286 catch (const AipsError &ae) {
1287 throw;
[881]1288 }
[368]1289 catch (const exception &ex) {
[881]1290 throw AipsError(String("STLineFinder::keepStrongestOnly - STL error: ")+ex.what());
[368]1291 }
1292
1293}
1294
[352]1295//
1296///////////////////////////////////////////////////////////////////////////////
1297
1298
1299///////////////////////////////////////////////////////////////////////////////
1300//
1301// LFLineListOperations - a class incapsulating operations with line lists
1302// The LF prefix stands for Line Finder
1303//
1304
[331]1305// concatenate two lists preserving the order. If two lines appear to
1306// be adjacent, they are joined into the new one
[352]1307void LFLineListOperations::addNewSearchResult(const std::list<pair<int, int> > &newlines,
[881]1308 std::list<std::pair<int, int> > &lines_list)
[331]1309 throw(AipsError)
1310{
1311 try {
1312 for (std::list<pair<int,int> >::const_iterator cli=newlines.begin();
1313 cli!=newlines.end();++cli) {
[881]1314
[996]1315 // the first item, which has a non-void intersection or touches
1316 // the new line
1317 std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(),
1318 lines_list.end(), IntersectsWith(*cli));
1319 // the last such item
1320 std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg,
1321 lines_list.end(), not1(IntersectsWith(*cli)));
[343]1322
1323 // extract all lines which intersect or touch a new one into
[996]1324 // a temporary buffer. This may invalidate the iterators
1325 // line_buffer may be empty, if no lines intersects with a new
1326 // one.
1327 std::list<pair<int,int> > lines_buffer;
1328 lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end);
[343]1329
[996]1330 // build a union of all intersecting lines
1331 pair<int,int> union_line=for_each(lines_buffer.begin(),
1332 lines_buffer.end(),BuildUnion(*cli)).result();
[881]1333
[996]1334 // search for a right place for the new line (union_line) and add
1335 std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(),
1336 lines_list.end(), LaterThan(union_line));
1337 lines_list.insert(pos2insert,union_line);
[331]1338 }
1339 }
1340 catch (const AipsError &ae) {
1341 throw;
[881]1342 }
[331]1343 catch (const exception &ex) {
[352]1344 throw AipsError(String("LFLineListOperations::addNewSearchResult - STL error: ")+ex.what());
[331]1345 }
1346}
[344]1347
1348// extend all line ranges to the point where a value stored in the
1349// specified vector changes (e.g. value-mean change its sign)
1350// This operation is necessary to include line wings, which are below
1351// the detection threshold. If lines becomes adjacent, they are
1352// merged together. Any masked channel stops the extension
[352]1353void LFLineListOperations::searchForWings(std::list<std::pair<int, int> > &newlines,
1354 const casa::Vector<casa::Int> &signs,
[996]1355 const casa::Vector<casa::Bool> &mask,
1356 const std::pair<int,int> &edge) throw(casa::AipsError)
[344]1357{
1358 try {
1359 for (std::list<pair<int,int> >::iterator li=newlines.begin();
1360 li!=newlines.end();++li) {
[996]1361 // update the left hand side
1362 for (int n=li->first-1;n>=edge.first;--n) {
1363 if (!mask[n]) break;
1364 if (signs[n]==signs[li->first] && signs[li->first])
1365 li->first=n;
1366 else break;
1367 }
1368 // update the right hand side
1369 for (int n=li->second;n<edge.second;++n) {
1370 if (!mask[n]) break;
1371 if (signs[n]==signs[li->second-1] && signs[li->second-1])
1372 li->second=n;
1373 else break;
1374 }
[344]1375 }
1376 // need to search for possible mergers.
1377 std::list<std::pair<int, int> > result_buffer;
1378 addNewSearchResult(newlines,result_buffer);
1379 newlines.clear();
1380 newlines.splice(newlines.end(),result_buffer);
1381 }
1382 catch (const AipsError &ae) {
1383 throw;
[881]1384 }
[344]1385 catch (const exception &ex) {
[352]1386 throw AipsError(String("LFLineListOperations::extendLines - STL error: ")+ex.what());
[344]1387 }
1388}
[352]1389
1390//
1391///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.