manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Added methods to Histogram and Distribution to get back the number of inserted samples and their total sum

This commit is contained in:
Paolo Cignoni 2014-05-13 10:54:42 +00:00
parent 5c7949d261
commit db0a706384
1 changed files with 244 additions and 304 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

548
vcg/math/histogram.h
View File

@ -8,7 +8,7 @@
* \ * * \ *
* All rights reserved. * * All rights reserved. *
* * * *
* This program is free software; you can redistribute it and/or modify * * This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by * * it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or * * the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. * * (at your option) any later version. *
@ -20,70 +20,6 @@
* for more details. * * for more details. *
* * * *
****************************************************************************/ ****************************************************************************/
/****************************************************************************
History
$Log: not supported by cvs2svn $
Revision 1.18 2008/03/05 11:21:49 cignoni
Heavily refactored the whole structure.
Some interfaces have been changed. Be careful.
Revision 1.17 2008/02/29 12:15:06 cignoni
added maxcount
Revision 1.16 2006/11/28 21:29:21 cignoni
Re added typedef Histogramf and Histogramd
Revision 1.15 2006/11/28 09:47:42 corsini
add documentation
fix typo
Revision 1.14 2006/05/04 00:09:53 cignoni
minor change: removed unused vars
Revision 1.13 2006/03/29 09:25:47 zifnab1974
extra includes necessary for compilation of meshlab on AMD 64 with gcc 3.4.5
Revision 1.12 2006/01/12 13:12:54 callieri
in FileWrite, added file closing after finishing
Revision 1.11 2005/09/16 11:51:23 cignoni
removed signed/unsigned warning
Revision 1.10 2005/06/17 00:54:55 cignoni
Corrected small bug in SetRange (H was resized to n instead of n+1)
Revision 1.9 2005/06/14 14:27:00 ganovelli
added include of algorithm
Revision 1.8 2005/06/10 14:59:39 cignoni
Added include assert.h and cast to ScalarType for a pow in SetRange() function.
Revision 1.7 2005/06/09 14:19:55 cignoni
Added typedef Histogramf and Histogramd
Revision 1.6 2005/06/07 09:37:33 ponchio
Added fabs() to variance, which can sometime be negative in case
of rounding errors (and sqrt chokes on it).
Revision 1.5 2005/06/07 07:44:08 cignoni
Added Percentile and removed small bug in Add
Revision 1.4 2005/04/04 10:48:35 cignoni
Added missing functions Avg, rms etc, now fully (almost) functional
Revision 1.3 2005/03/14 09:23:40 cignoni
Added missing include<vector>
Revision 1.2 2004/08/25 15:15:26 ganovelli
minor changes to comply gcc compiler (typename's and stuff)
Revision 1.1 2004/06/24 09:12:28 cignoni
Initial Release
****************************************************************************/
#ifndef __VCG_HISTOGRAM #ifndef __VCG_HISTOGRAM
#define __VCG_HISTOGRAM #define __VCG_HISTOGRAM
@ -97,83 +33,86 @@ Initial Release
namespace vcg { namespace vcg {
template <class ScalarType> template <class ScalarType>
class Distribution class Distribution
{ {
private: private:
std::vector<ScalarType> vec; std::vector<ScalarType> vec;
bool dirty; bool dirty;
double valSum; double valSum;
double sqrdValSum; double sqrdValSum;
double avg; double avg;
double sqrdAvg; double sqrdAvg;
double rms; double rms;
double min_v; double min_v;
double max_v; double max_v;
public: public:
Distribution() { Clear(); } Distribution() { Clear(); }
void Clear() void Clear()
{ {
vec.clear(); vec.clear();
dirty=true; dirty=true;
min_v = std::numeric_limits<float>::max(); min_v = std::numeric_limits<float>::max();
max_v = -std::numeric_limits<float>::max(); max_v = -std::numeric_limits<float>::max();
} }
void Add(const ScalarType v)
{
vec.push_back(v);
dirty=true;
if(v<min_v) min_v=v;
if(v>max_v) max_v=v;
}
ScalarType Min() { return min_v; }
ScalarType Max() { return max_v; }
ScalarType Avg(){ DirtyCheck(); return avg;} void Add(const ScalarType v)
//! Returns the Root Mean Square of the data. {
ScalarType RMS(){ DirtyCheck(); return rms;} vec.push_back(v);
dirty=true;
//! Returns the variance of the data. if(v<min_v) min_v=v;
// the average of the squares less the square of the average. if(v>max_v) max_v=v;
ScalarType Variance(){ DirtyCheck(); return sqrdAvg - avg*avg ;} }
//! Returns the standard deviation of the data. ScalarType Min() const { return min_v; }
ScalarType StandardDeviation(){ DirtyCheck(); return sqrt( Variance() );} ScalarType Max() const { return max_v; }
void DirtyCheck() ScalarType Cnt() const { return ScalarType(vec.size()); }
{
if(!dirty) return; ScalarType Sum(){ DirtyCheck(); return valSum; }
std::sort(vec.begin(),vec.end()); ScalarType Avg(){ DirtyCheck(); return avg;}
valSum=0; //! Returns the Root Mean Square of the data.
sqrdValSum=0; ScalarType RMS(){ DirtyCheck(); return rms;}
typename std::vector<ScalarType>::iterator vi;
for(vi=vec.begin();vi!=vec.end();++vi) //! \brief Returns the variance of the data.
{ /// the average of the squares less the square of the average.
valSum += double(*vi); ScalarType Variance(){ DirtyCheck(); return sqrdAvg - avg*avg ;}
sqrdValSum += double(*vi)*double(*vi);
} //! Returns the standard deviation of the data.
avg = valSum/double(vec.size()); ScalarType StandardDeviation(){ DirtyCheck(); return sqrt( Variance() );}
sqrdAvg = sqrdValSum/double(vec.size());
rms = math::Sqrt(sqrdAvg); void DirtyCheck()
dirty=false; {
} if(!dirty) return;
std::sort(vec.begin(),vec.end());
ScalarType Percentile(ScalarType perc) valSum=0;
{ sqrdValSum=0;
assert(!vec.empty()); typename std::vector<ScalarType>::iterator vi;
assert(perc>=0 && perc<=1); for(vi=vec.begin();vi!=vec.end();++vi)
DirtyCheck(); {
int index = vec.size() *perc -1; valSum += double(*vi);
sqrdValSum += double(*vi)*double(*vi);
if(index< 0 ) index = 0; }
avg = valSum/double(vec.size());
return vec[index]; sqrdAvg = sqrdValSum/double(vec.size());
} rms = math::Sqrt(sqrdAvg);
dirty=false;
}
ScalarType Percentile(ScalarType perc)
{
assert(!vec.empty());
assert(perc>=0 && perc<=1);
DirtyCheck();
int index = vec.size() *perc -1;
if(index< 0 ) index = 0;
return vec[index];
}
}; };
@ -183,17 +122,17 @@ public:
* *
* This class implements a single-value histogram. * This class implements a single-value histogram.
*/ */
template <class ScalarType> template <class ScalarType>
class Histogram class Histogram
{ {
// public data members // public data members
protected: protected:
std::vector <ScalarType> H; //! Counters for bins. std::vector <ScalarType> H; //! Counters for bins.
std::vector <ScalarType> R; //! Range for bins. std::vector <ScalarType> R; //! Range for bins.
ScalarType minv; //! Minimum value. ScalarType minv; //! Minimum value.
ScalarType maxv; //! Maximum value. ScalarType maxv; //! Maximum value.
ScalarType minElem; //! Minimum value. ScalarType minElem; //! Minimum value.
ScalarType maxElem; //! Maximum value. ScalarType maxElem; //! Maximum value.
int n; //! Number of vaild intervals stored between minv and maxv. int n; //! Number of vaild intervals stored between minv and maxv.
@ -201,91 +140,92 @@ protected:
/// incrementally updated values /// incrementally updated values
ScalarType cnt; //! Number of accumulated samples. ScalarType cnt; //! Number of accumulated samples.
ScalarType avg; //! Average. ScalarType sum; //! Average.
ScalarType rms; //! Root mean square. ScalarType rms; //! Root mean square.
/** /**
* Returns the index of the bin which contains a given value. * Returns the index of the bin which contains a given value.
*/ */
int BinIndex(ScalarType val) ; int BinIndex(ScalarType val) ;
// public methods // public methods
public: public:
/**
* Set the histogram values.
*
* This method is used to correctly initialize the bins of the histogram.
* n is the number of valid intervals between minv and maxv.
* for a more robust working, the Histogram class stores also the two out of range intervals (-inf, minv] and [maxv, +inf)
* Each bin is left closed (eg it contains the value
* The \a gamma parameter is applied to modify the distribution of the ranges of the bins. Default uniform distibution.
*
*/
void SetRange(ScalarType _minv, ScalarType _maxv, int _n,ScalarType gamma=1.0 );
ScalarType MinV() {return minv;}; //! Minimum value. /**
ScalarType MaxV() {return maxv;}; //! Maximum value. * Set the histogram values.
*
* This method is used to correctly initialize the bins of the histogram.
* n is the number of valid intervals between minv and maxv.
* for a more robust working, the Histogram class stores also the two out of range intervals (-inf, minv] and [maxv, +inf)
* Each bin is left closed (eg it contains the value
* The \a gamma parameter is applied to modify the distribution of the ranges of the bins. Default uniform distibution.
*
*/
void SetRange(ScalarType _minv, ScalarType _maxv, int _n,ScalarType gamma=1.0 );
ScalarType MinV() {return minv;} //! Minimum value.
ScalarType MinElem() {return minElem;}; //! Minimum element added to the histogram. It could be < or > than MinV;. ScalarType MaxV() {return maxv;} //! Maximum value.
ScalarType MaxElem() {return maxElem;}; //! Maximum element added to the histogram. It could be < or > than MinV;.. ScalarType Sum() {return sum;} //! Total sum of inserted values.
ScalarType Cnt() {return cnt;}
/** ScalarType MinElem() {return minElem;} //! Minimum element added to the histogram. It could be < or > than MinV;.
* Add a new value to the histogram. ScalarType MaxElem() {return maxElem;} //! Maximum element added to the histogram. It could be < or > than MinV;..
*
* The statistics related to the histogram data (average, RMS, etc.) are /**
* also updated. * Add a new value to the histogram.
*/ *
* The statistics related to the histogram data (average, RMS, etc.) are
* also updated.
*/
void Add(ScalarType v, ScalarType increment=ScalarType(1.0)); void Add(ScalarType v, ScalarType increment=ScalarType(1.0));
ScalarType MaxCount() const; ScalarType MaxCount() const;
int BinNum() const {return n;}; int BinNum() const {return n;}
ScalarType BinCount(ScalarType v); ScalarType BinCount(ScalarType v);
ScalarType BinCountInd(int index) {return H[index];} ScalarType BinCountInd(int index) {return H[index];}
ScalarType BinCount(ScalarType v, ScalarType width); ScalarType BinCount(ScalarType v, ScalarType width);
ScalarType BinLowerBound(int index) {return R[index];} ScalarType BinLowerBound(int index) {return R[index];}
ScalarType BinUpperBound(int index) {return R[index+1];}; ScalarType BinUpperBound(int index) {return R[index+1];}
ScalarType RangeCount(ScalarType rangeMin, ScalarType rangeMax); ScalarType RangeCount(ScalarType rangeMin, ScalarType rangeMax);
ScalarType BinWidth(ScalarType v); ScalarType BinWidth(ScalarType v);
/** /**
* Returns the value corresponding to a given percentile of the data. * Returns the value corresponding to a given percentile of the data.
* *
* The percentile range between 0 and 1. * The percentile range between 0 and 1.
*/ */
ScalarType Percentile(ScalarType frac) const; ScalarType Percentile(ScalarType frac) const;
//! Returns the average of the data.
ScalarType Avg(){ return avg/cnt;}
//! Returns the Root Mean Square of the data.
ScalarType RMS(){ return sqrt(rms/double(cnt));}
//! Returns the variance of the data.
ScalarType Variance(){ return fabs(rms/cnt-Avg()*Avg());}
//! Returns the standard deviation of the data.
ScalarType StandardDeviation(){ return sqrt(Variance());}
//! Dump the histogram to a file. //! Returns the average of the data.
void FileWrite(const std::string &filename); ScalarType Avg(){ return sum/cnt;}
//! Reset histogram data. //! Returns the Root Mean Square of the data.
void Clear(); ScalarType RMS(){ return sqrt(rms/double(cnt));}
//! Returns the variance of the data.
ScalarType Variance(){ return fabs(rms/cnt-Avg()*Avg());}
//! Returns the standard deviation of the data.
ScalarType StandardDeviation(){ return sqrt(Variance());}
//! Dump the histogram to a file.
void FileWrite(const std::string &filename);
//! Reset histogram data.
void Clear();
}; };
template <class ScalarType> template <class ScalarType>
void Histogram<ScalarType>::Clear() void Histogram<ScalarType>::Clear()
{ {
H.clear(); H.clear();
R.clear(); R.clear();
cnt=0; cnt=0;
avg=0; sum=0;
rms=0; rms=0;
n=0; n=0;
minv=0; minv=0;
maxv=1; maxv=1;
minElem = std::numeric_limits<ScalarType>::max(); minElem = std::numeric_limits<ScalarType>::max();
maxElem = -std::numeric_limits<ScalarType>::max(); maxElem = -std::numeric_limits<ScalarType>::max();
} }
@ -293,169 +233,169 @@ void Histogram<ScalarType>::Clear()
/* /*
Note that the histogram holds <n> valid bins plus two semi-infinite bins. Note that the histogram holds <n> valid bins plus two semi-infinite bins.
R[0] = -inf R[0] = -inf
R[1] = minv R[1] = minv
R[n+1] = maxv R[n+1] = maxv
R[n+2] = +inf R[n+2] = +inf
Eg. SetRange(0, 10, 5) asks for 5 intervals covering the 0..10 range Eg. SetRange(0, 10, 5) asks for 5 intervals covering the 0..10 range
H[0] H[1] H[2] H[3] H[4] H[5] H[6] H[0] H[1] H[2] H[3] H[4] H[5] H[6]
-inf 0 2 4 6 8 10 +inf -inf 0 2 4 6 8 10 +inf
R[0] R[1] R[2] R[3] R[4] R[5] R[6] R[7] R[0] R[1] R[2] R[3] R[4] R[5] R[6] R[7]
*/ */
template <class ScalarType> template <class ScalarType>
void Histogram<ScalarType>::SetRange(ScalarType _minv, ScalarType _maxv, int _n, ScalarType gamma) void Histogram<ScalarType>::SetRange(ScalarType _minv, ScalarType _maxv, int _n, ScalarType gamma)
{ {
// reset data // reset data
Clear(); Clear();
minv=_minv;maxv=_maxv;n=_n;
H.resize(n+2);
fill(H.begin(),H.end(),0);
R.resize(n+3);
R[0] = - std::numeric_limits< ScalarType >::max();
R[n+2] = std::numeric_limits< ScalarType >::max();
double delta=(maxv-minv); minv=_minv;maxv=_maxv;n=_n;
if(gamma==1) H.resize(n+2);
{ fill(H.begin(),H.end(),0);
for(int i=0; i<=n; ++i) R.resize(n+3);
R[i+1] = minv + delta*ScalarType(i)/n;
} R[0] = - std::numeric_limits< ScalarType >::max();
else R[n+2] = std::numeric_limits< ScalarType >::max();
{
for(int i=0; i<=n; ++i) double delta=(maxv-minv);
R[i+1] = minv + delta*pow(ScalarType(i)/n,gamma); if(gamma==1)
} {
for(int i=0; i<=n; ++i)
R[i+1] = minv + delta*ScalarType(i)/n;
}
else
{
for(int i=0; i<=n; ++i)
R[i+1] = minv + delta*pow(ScalarType(i)/n,gamma);
}
} }
template <class ScalarType> template <class ScalarType>
int Histogram<ScalarType>::BinIndex(ScalarType val) int Histogram<ScalarType>::BinIndex(ScalarType val)
{ {
// lower_bound returns the furthermost iterator i in [first, last) such that, for every iterator j in [first, i), *j < value. // lower_bound returns the furthermost iterator i in [first, last) such that, for every iterator j in [first, i), *j < value.
// E.g. An iterator pointing to the first element "not less than" val, or end() if every element is less than val. // E.g. An iterator pointing to the first element "not less than" val, or end() if every element is less than val.
typename std::vector<ScalarType>::iterator it = lower_bound(R.begin(),R.end(),val); typename std::vector<ScalarType>::iterator it = lower_bound(R.begin(),R.end(),val);
assert(it!=R.begin()); assert(it!=R.begin());
assert(it!=R.end()); assert(it!=R.end());
assert((*it)>=val); assert((*it)>=val);
int pos = it-R.begin(); int pos = it-R.begin();
assert(pos >=1); assert(pos >=1);
pos -= 1; pos -= 1;
assert (R[pos] < val); assert (R[pos] < val);
assert ( val <= R[pos+1] ); assert ( val <= R[pos+1] );
return pos; return pos;
} }
/* /*
H[0] H[1] H[2] H[3] H[4] H[5] H[6] H[0] H[1] H[2] H[3] H[4] H[5] H[6]
-inf 0 2 4 6 8 10 +inf -inf 0 2 4 6 8 10 +inf
R[0] R[1] R[2] R[3] R[4] R[5] R[6] R[7] R[0] R[1] R[2] R[3] R[4] R[5] R[6] R[7]
asking for 3.14 lower bound will return an iterator pointing to R[3]==4; and will increase H[2] asking for 3.14 lower bound will return an iterator pointing to R[3]==4; and will increase H[2]
asking for 4 lower bound will return an iterator pointing to R[3]==4; and will increase H[2] asking for 4 lower bound will return an iterator pointing to R[3]==4; and will increase H[2]
*/ */
template <class ScalarType> template <class ScalarType>
void Histogram<ScalarType>::Add(ScalarType v, ScalarType increment) void Histogram<ScalarType>::Add(ScalarType v, ScalarType increment)
{ {
int pos=BinIndex(v); int pos=BinIndex(v);
if(v<minElem) minElem=v; if(v<minElem) minElem=v;
if(v>maxElem) maxElem=v; if(v>maxElem) maxElem=v;
assert((pos>=0)&&(pos<=n+1)); assert((pos>=0)&&(pos<=n+1));
H[pos]+=increment; H[pos]+=increment;
cnt+=increment; cnt+=increment;
avg+=v*increment; sum+=v*increment;
rms += (v*v)*increment; rms += (v*v)*increment;
} }
template <class ScalarType> template <class ScalarType>
ScalarType Histogram<ScalarType>::BinCount(ScalarType v) ScalarType Histogram<ScalarType>::BinCount(ScalarType v)
{ {
return H[BinIndex(v)]; return H[BinIndex(v)];
} }
template <class ScalarType> template <class ScalarType>
ScalarType Histogram<ScalarType>::BinCount(ScalarType v, ScalarType width) ScalarType Histogram<ScalarType>::BinCount(ScalarType v, ScalarType width)
{ {
return RangeCount(v-width/2.0,v+width/2.0); return RangeCount(v-width/2.0,v+width/2.0);
} }
template <class ScalarType> template <class ScalarType>
ScalarType Histogram<ScalarType>::RangeCount(ScalarType rangeMin, ScalarType rangeMax) ScalarType Histogram<ScalarType>::RangeCount(ScalarType rangeMin, ScalarType rangeMax)
{ {
int firstBin=BinIndex(rangeMin); int firstBin=BinIndex(rangeMin);
int lastBin=BinIndex (rangeMax); int lastBin=BinIndex (rangeMax);
ScalarType sum=0; ScalarType sum=0;
for(int i=firstBin; i<=lastBin;++i) for(int i=firstBin; i<=lastBin;++i)
sum+=H[i]; sum+=H[i];
return sum; return sum;
} }
template <class ScalarType> template <class ScalarType>
ScalarType Histogram<ScalarType>::BinWidth(ScalarType v) ScalarType Histogram<ScalarType>::BinWidth(ScalarType v)
{ {
int pos=BinIndex(v); int pos=BinIndex(v);
return R[pos+1]-R[pos]; return R[pos+1]-R[pos];
} }
template <class ScalarType> template <class ScalarType>
void Histogram<ScalarType>::FileWrite(const std::string &filename) void Histogram<ScalarType>::FileWrite(const std::string &filename)
{ {
FILE *fp; FILE *fp;
fp=fopen(filename.c_str(),"w"); fp=fopen(filename.c_str(),"w");
for(unsigned int i=0; i<H.size(); i++)
fprintf (fp,"%12.8lf , %12.8lf \n",R[i],double(H[i])/cnt);
fclose(fp); for(unsigned int i=0; i<H.size(); i++)
fprintf (fp,"%12.8lf , %12.8lf \n",R[i],double(H[i])/cnt);
fclose(fp);
} }
template <class ScalarType> template <class ScalarType>
ScalarType Histogram<ScalarType>::MaxCount() const ScalarType Histogram<ScalarType>::MaxCount() const
{ {
return *(std::max_element(H.begin(),H.end())); return *(std::max_element(H.begin(),H.end()));
} }
// Return the scalar value <r> such that there are <frac> samples <= <r>. // Return the scalar value <r> such that there are <frac> samples <= <r>.
// E.g. Percentile(0.0) will return R[1] e.g. min value // E.g. Percentile(0.0) will return R[1] e.g. min value
// E.g. Percentile(1.0) will return R[n+1] e.g max value // E.g. Percentile(1.0) will return R[n+1] e.g max value
template <class ScalarType> template <class ScalarType>
ScalarType Histogram<ScalarType>::Percentile(ScalarType frac) const ScalarType Histogram<ScalarType>::Percentile(ScalarType frac) const
{ {
if(H.size()==0 && R.size()==0) if(H.size()==0 && R.size()==0)
return 0; return 0;
// check percentile range
assert(frac >= 0 && frac <= 1);
ScalarType sum=0,partsum=0;
size_t i;
// useless summation just to be sure // check percentile range
for(i=0;i<H.size();i++) sum+=H[i]; assert(frac >= 0 && frac <= 1);
assert(sum==cnt);
sum*=frac; ScalarType sum=0,partsum=0;
for(i=0; i<H.size(); i++) size_t i;
{
partsum+=H[i]; // useless summation just to be sure
if(partsum>=sum) break; for(i=0;i<H.size();i++) sum+=H[i];
} assert(sum==cnt);
assert(i<H.size()); sum*=frac;
for(i=0; i<H.size(); i++)
return R[i+1]; {
partsum+=H[i];
if(partsum>=sum) break;
}
assert(i<H.size());
return R[i+1];
} }
typedef Histogram<double> Histogramd ; typedef Histogram<double> Histogramd ;