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