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