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