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