Changeset 1642 for trunk


Ignore:
Timestamp:
10/03/09 00:13:24 (15 years ago)
Author:
Max Voronkov
Message:

added a new helper class to the line finder (compilable, but not yet used). It will allow to improve line finder by adding more options of noise estimation

Location:
trunk/src
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/STLineFinder.cpp

    r1641 r1642  
    3434#include "STLineFinder.h"
    3535#include "STFitter.h"
     36#include "IndexedCompare.h"
    3637
    3738// STL
     
    200201///////////////////////////////////////////////////////////////////////////////
    201202
     203///////////////////////////////////////////////////////////////////////////////
     204//
     205// LFNoiseEstimator  a helper class designed to estimate off-line variance
     206//                   using statistics depending on the distribution of
     207//                   values (e.g. like a median)
     208//
     209//                   Two statistics are supported: median and an average of
     210//                   80% of smallest values.
     211//
     212
     213struct LFNoiseEstimator {
     214   // construct an object
     215   // size - maximum sample size. After a size number of elements is processed
     216   // any new samples would cause the algorithm to drop the oldest samples in the
     217   // buffer.
     218   explicit LFNoiseEstimator(size_t size);
     219
     220   // add a new sample
     221   // in - the new value
     222   void add(float in);
     223
     224   // median of the distribution
     225   float median() const;
     226
     227   // mean of lowest 80% of the samples
     228   float meanLowest80Percent() const;
     229
     230protected:
     231   // update cache of sorted indices
     232   // (it is assumed that itsSampleNumber points to the newly
     233   // replaced element)
     234   void updateSortedCache() const;
     235
     236   // build sorted cache from the scratch
     237   void buildSortedCache() const;
     238
     239   // number of samples accumulated so far
     240   // (can be less than the buffer size)
     241   size_t numberOfSamples() const;
     242
     243   // this helper method builds the cache if
     244   // necessary using one of the methods
     245   void fillCacheIfNecessary() const;
     246
     247private:
     248   // buffer with samples (unsorted)
     249   std::vector<float> itsVariances;
     250   // current sample number (<=itsVariances.size())
     251   size_t itsSampleNumber;
     252   // true, if the buffer all values in the sample buffer are used
     253   bool itsBufferFull;
     254   // cached indices into vector of samples
     255   mutable std::vector<size_t> itsSortedIndices;
     256   // true if any of the statistics have been obtained at least
     257   // once. This flag allows to implement a more efficient way of
     258   // calculating statistics, if they are needed at once and not
     259   // after each addition of a new element
     260   mutable bool itsStatisticsAccessed;
     261};
     262
     263//
     264///////////////////////////////////////////////////////////////////////////////
     265
     266
    202267} // namespace asap
     268
     269///////////////////////////////////////////////////////////////////////////////
     270//
     271// LFNoiseEstimator  a helper class designed to estimate off-line variance
     272//                   using statistics depending on the distribution of
     273//                   values (e.g. like a median)
     274//
     275//                   Two statistics are supported: median and an average of
     276//                   80% of smallest values.
     277//
     278
     279// construct an object
     280// size - maximum sample size. After a size number of elements is processed
     281// any new samples would cause the algorithm to drop the oldest samples in the
     282// buffer.
     283LFNoiseEstimator::LFNoiseEstimator(size_t size) : itsVariances(size),
     284     itsSampleNumber(0), itsBufferFull(false), itsSortedIndices(size),
     285     itsStatisticsAccessed(false)
     286{
     287   AlwaysAssert(size>0,AipsError);
     288}
     289
     290
     291// add a new sample
     292// in - the new value
     293void LFNoiseEstimator::add(float in)
     294{
     295   itsVariances[itsSampleNumber] = in;
     296
     297   if (itsStatisticsAccessed) {
     298       // only do element by element addition if on-the-fly
     299       // statistics are needed
     300       updateSortedCache();
     301   }
     302
     303   // advance itsSampleNumber now
     304   ++itsSampleNumber;
     305   if (itsSampleNumber == itsVariances.size()) {
     306       itsSampleNumber = 0;
     307       itsBufferFull = true;
     308   }
     309   AlwaysAssert(itsSampleNumber<itsVariances.size(),AipsError);
     310}
     311
     312// number of samples accumulated so far
     313// (can be less than the buffer size)
     314size_t LFNoiseEstimator::numberOfSamples() const
     315{
     316  // the number of samples accumulated so far may be less than the
     317  // buffer size
     318  const size_t nSamples = itsBufferFull ? itsVariances.size(): itsSampleNumber;
     319  AlwaysAssert( (nSamples > 0) && (nSamples < itsVariances.size()), AipsError);
     320  return nSamples;
     321}
     322
     323// this helper method builds the cache if
     324// necessary using one of the methods
     325void LFNoiseEstimator::fillCacheIfNecessary() const
     326{
     327  if (!itsStatisticsAccessed) {
     328      if ((itsSampleNumber!=0) || itsBufferFull) {
     329          // build the whole cache efficiently
     330          buildSortedCache();
     331      } else {
     332          updateSortedCache();
     333      }
     334      itsStatisticsAccessed = true;
     335  } // otherwise, it is updated in 'add' using on-the-fly method
     336}
     337
     338// median of the distribution
     339float LFNoiseEstimator::median() const
     340{
     341  fillCacheIfNecessary();
     342  // the number of samples accumulated so far may be less than the
     343  // buffer size
     344  const size_t nSamples = numberOfSamples();
     345  const size_t medSample = nSamples / 2;
     346  AlwaysAssert(medSample < itsSortedIndices.size(), AipsError);
     347  return itsVariances[itsSortedIndices[medSample]];
     348}
     349
     350// mean of lowest 80% of the samples
     351float LFNoiseEstimator::meanLowest80Percent() const
     352{
     353  fillCacheIfNecessary();
     354  // the number of samples accumulated so far may be less than the
     355  // buffer size
     356  const size_t nSamples = numberOfSamples();
     357  float result = 0;
     358  size_t numpt=size_t(0.8*nSamples);
     359  if (!numpt) {
     360      numpt=nSamples; // no much else left,
     361                     // although it is very inaccurate
     362  }
     363  AlwaysAssert( (numpt > 0) && (numpt<itsSortedIndices.size()), AipsError);
     364  for (size_t ch=0; ch<numpt; ++ch) {
     365       result += itsVariances[itsSortedIndices[ch]];
     366  }
     367  result /= float(numpt);
     368  return result;
     369}
     370
     371// update cache of sorted indices
     372// (it is assumed that itsSampleNumber points to the newly
     373// replaced element)
     374void LFNoiseEstimator::updateSortedCache() const
     375{
     376  // the number of samples accumulated so far may be less than the
     377  // buffer size
     378  const size_t nSamples = numberOfSamples();
     379
     380  if (itsBufferFull) {
     381      // first find the index of the element which is being replaced
     382      size_t index = nSamples;
     383      for (size_t i=0; i<nSamples; ++i) {
     384           AlwaysAssert(i < itsSortedIndices.size(), AipsError);
     385           if (itsSortedIndices[i] == itsSampleNumber) {
     386               index = i;
     387               break;
     388           }
     389      }
     390      AlwaysAssert( index < nSamples, AipsError);
     391
     392      const vector<size_t>::iterator indStart = itsSortedIndices.begin();
     393      // merge this element with preceeding block first
     394      if (index != 0) {
     395          // merge indices on the basis of variances
     396          inplace_merge(indStart,indStart+index,indStart+index+1,
     397                        indexedCompare<size_t>(itsVariances.begin()));
     398      }
     399      // merge with the following block
     400      if (index + 1 != nSamples) {
     401          // merge indices on the basis of variances
     402          inplace_merge(indStart,indStart+index+1,indStart+nSamples,
     403                        indexedCompare<size_t>(itsVariances.begin()));
     404      }
     405  } else {
     406      // itsSampleNumber is the index of the new element
     407      AlwaysAssert(itsSampleNumber < itsSortedIndices.size(), AipsError);
     408      itsSortedIndices[itsSampleNumber] = itsSampleNumber;
     409      if (itsSampleNumber >= 1) {
     410          // we have to place this new sample in
     411          const vector<size_t>::iterator indStart = itsSortedIndices.begin();
     412          // merge indices on the basis of variances
     413          inplace_merge(indStart,indStart+itsSampleNumber,indStart+itsSampleNumber+1,
     414                        indexedCompare<size_t>(itsVariances.begin()));
     415      }
     416  }
     417}
     418
     419// build sorted cache from the scratch
     420void LFNoiseEstimator::buildSortedCache() const
     421{
     422  // the number of samples accumulated so far may be less than the
     423  // buffer size
     424  const size_t nSamples = numberOfSamples();
     425  AlwaysAssert(nSamples < itsSortedIndices.size(), AipsError);
     426  for (size_t i=0; i<nSamples; ++i) {
     427       itsSortedIndices[i]=i;
     428  }
     429
     430  // sort indices, but check the array of variances
     431  const vector<size_t>::iterator indStart = itsSortedIndices.begin();
     432  stable_sort(indStart,indStart+nSamples, indexedCompare<size_t>(itsVariances.begin()));
     433}
     434
     435//
     436///////////////////////////////////////////////////////////////////////////////
    203437
    204438///////////////////////////////////////////////////////////////////////////////
Note: See TracChangeset for help on using the changeset viewer.