histogram.hh
Go to the documentation of this file.
/* -*- mia-c++ -*-
*
* This file is part of MIA - a toolbox for medical image analysis
* Copyright (c) Leipzig, Madrid 1999-2017 Gert Wollny
*
* MIA is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with MIA; if not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef mia_core_histogram_hh
#define mia_core_histogram_hh
#include <mia/core/defines.hh>
#include <cmath>
#include <cassert>
#include <vector>
#include <mia/core/defines.hh>
#include <mia/core/msgstream.hh>
#include <boost/type_traits.hpp>
NS_MIA_BEGIN
template <typename T>
class THistogramFeeder
{
public:
typedef T value_type;
THistogramFeeder(T min, T max, size_t bins);
size_t size() const;
size_t index(T x) const;
T value(size_t k) const;
private:
T m_min;
T m_max;
size_t m_bins;
double m_step;
double m_inv_step;
};
template <>
class THistogramFeeder<unsigned char>
{
public:
typedef unsigned char value_type;
THistogramFeeder(unsigned char min, unsigned char max, size_t bins);
size_t size() const;
size_t index(unsigned char x) const;
unsigned char value(size_t k) const;
};
typedef THistogramFeeder<unsigned char> CUBHistogramFeeder;
template <typename Feeder>
class THistogram
{
public:
typedef std::vector<size_t>::const_iterator const_iterator;
typedef std::pair<typename Feeder::value_type, size_t> value_type;
typedef std::pair<typename Feeder::value_type, typename Feeder::value_type> range_type;
THistogram(const Feeder& f);
THistogram(const THistogram<Feeder>& org, double perc);
void push(typename Feeder::value_type x);
void push(typename Feeder::value_type x, size_t count);
template <typename Iterator>
void push_range(Iterator begin, Iterator end);
size_t size() const;
const_iterator begin() const;
const_iterator end() const;
size_t operator [] (size_t idx) const;
const value_type at(size_t idx) const;
typename Feeder::value_type median() const;
typename Feeder::value_type MAD() const;
double average() const;
double deviation() const;
double excess_kurtosis() const;
double skewness() const;
range_type get_reduced_range(double remove) const;
private:
Feeder m_feeder;
std::vector<size_t> m_histogram;
size_t m_n;
};
// inline inplementation
template <typename T>
THistogramFeeder<T>::THistogramFeeder(T min, T max, size_t bins):
m_min(min),
m_max(max),
m_bins(bins),
m_step(( double(max) - double(min) ) / double(bins - 1)),
m_inv_step(double(bins - 1) / (double(max) - double(min)))
{
}
template <typename T>
size_t THistogramFeeder<T>::size() const
{
return m_bins;
}
template <typename T>
inline size_t THistogramFeeder<T>::index(T x) const
{
double val = floor(m_inv_step * (x - m_min) + 0.5);
if (val < 0)
return 0;
if (val < m_bins)
return val;
return m_bins - 1;
}
template <typename T>
T THistogramFeeder<T>::value(size_t k) const
{
return k * m_step + m_min;
}
inline THistogramFeeder<unsigned char>::THistogramFeeder(unsigned char /*min*/, unsigned char /*max*/, size_t /*bins*/)
{
}
inline size_t THistogramFeeder<unsigned char>::size() const
{
return 256;
}
inline
size_t THistogramFeeder<unsigned char>::index(unsigned char x) const
{
return x;
}
inline
unsigned char THistogramFeeder<unsigned char>::value(size_t k) const
{
return k;
}
template <typename Feeder>
THistogram<Feeder>::THistogram(const Feeder& f):
m_feeder(f),
m_histogram(f.size()),
m_n(0)
{
}
template <typename Feeder>
THistogram<Feeder>::THistogram(const THistogram<Feeder>& org, double perc):
m_feeder(org.m_feeder),
m_histogram(m_feeder.size()),
m_n(0)
{
size_t n = (size_t)(org.m_n * (1.0 - perc));
size_t i = 0;
while (n > m_n && i < m_histogram.size()) {
m_n += org.m_histogram[i];
m_histogram[i] = org.m_histogram[i];
++i;
}
}
template <typename Feeder>
size_t THistogram<Feeder>::size() const
{
return m_histogram.size();
}
template <typename Feeder>
void THistogram<Feeder>::push(typename Feeder::value_type x)
{
++m_n;
++m_histogram[m_feeder.index(x)];
}
template <typename Feeder>
template <typename Iterator>
void THistogram<Feeder>::push_range(Iterator begin, Iterator end)
{
while (begin != end)
push(*begin++);
}
template <typename Feeder>
void THistogram<Feeder>::push(typename Feeder::value_type x, size_t count)
{
m_n += count;
m_histogram[m_feeder.index(x)] += count;
}
template <typename Feeder>
typename THistogram<Feeder>::const_iterator THistogram<Feeder>::begin() const
{
return m_histogram.begin();
}
template <typename Feeder>
typename THistogram<Feeder>::const_iterator THistogram<Feeder>::end() const
{
return m_histogram.end();
}
template <typename Feeder>
size_t THistogram<Feeder>::operator [] (size_t idx) const
{
assert(idx < m_histogram.size());
return m_histogram[idx];
}
template <typename Feeder>
typename Feeder::value_type THistogram<Feeder>::median() const
{
float n_2 = m_n / 2.0f;
float sum = 0;
size_t k = 0;
while ( sum < n_2 )
sum += m_histogram[k++];
return m_feeder.value(k > 0 ? k - 1 : k);
}
template <typename Feeder>
typename Feeder::value_type THistogram<Feeder>::MAD() const
{
typedef typename Feeder::value_type T;
T m = median();
THistogram<Feeder> help(m_feeder);
;
for (size_t k = 0; k < size(); ++k) {
T v = m_feeder.value(k);
help.push(v > m ? v - m : m - v, m_histogram[k]);
}
return help.median();
}
template <typename Feeder>
const typename THistogram<Feeder>::value_type THistogram<Feeder>::at(size_t idx) const
{
if (idx < m_histogram.size())
return value_type(m_feeder.value(idx), m_histogram[idx]);
else
return value_type(m_feeder.value(idx), 0);
}
template <typename Feeder>
double THistogram<Feeder>::average() const
{
if (m_n < 1)
return 0.0;
double sum = 0.0;
for (size_t i = 0; i < size(); ++i) {
const typename THistogram<Feeder>::value_type value = at(i);
sum += value.first * value.second;
}
return sum / m_n;
}
template <typename Feeder>
double THistogram<Feeder>::excess_kurtosis() const
{
double mu = average();
double sum1 = 0.0;
double sum2 = 0.0;
double sum3 = 0.0;
for (size_t i = 0; i < size(); ++i) {
const auto value = at(i);
sum1 += value.first * value.second;
sum2 += value.first * value.first * value.second;
double h = (value.first - mu);
h *= h;
h *= h;
sum3 += h * value.second;
}
double sigma2 = (sum2 - sum1 * sum1 / m_n) / m_n;
double mu4 = sum3 / m_n;
return mu4 / (sigma2 * sigma2) - 3;
}
template <typename Feeder>
double THistogram<Feeder>::skewness() const
{
double mu = average();
double sum1 = 0.0;
double sum2 = 0.0;
double sum3 = 0.0;
for (size_t i = 0; i < size(); ++i) {
const auto value = at(i);
sum1 += value.first * value.second;
sum2 += value.first * value.first * value.second;
double h = (value.first - mu);
sum3 += h * h * h * value.second;
}
double sigma2 = (sum2 - sum1 * sum1 / m_n) / m_n;
double mu4 = sum3 / m_n;
return mu4 / (sigma2 * sqrt(sigma2));
}
template <typename Feeder>
double THistogram<Feeder>::deviation() const
{
if (m_n < 2)
return 0.0;
double sum = 0.0;
double sum2 = 0.0;
for (size_t i = 0; i < size(); ++i) {
const typename THistogram<Feeder>::value_type value = at(i);
sum += value.first * value.second;
sum2 += value.first * value.first * value.second;
}
return sqrt((sum2 - sum * sum / m_n) / (m_n - 1));
}
template <typename Feeder>
typename THistogram<Feeder>::range_type
THistogram<Feeder>::get_reduced_range(double remove) const
{
assert(remove >= 0.0 && remove < 49.0);
long remove_count = static_cast<long>(remove * m_n / 100.0);
range_type result(m_feeder.value(0), m_feeder.value(m_histogram.size() - 1));
if (remove_count > 0) {
long low_end = -1;
long counted_pixels_low = 0;
while (counted_pixels_low < remove_count && low_end < (long)m_histogram.size())
counted_pixels_low += m_histogram[++low_end];
result.first = m_feeder.value(low_end);
long high_end = m_histogram.size();
long counted_pixels_high = 0;
while (counted_pixels_high <= remove_count && high_end > 0)
counted_pixels_high += m_histogram[--high_end];
cvdebug() << " int range = " << low_end << ", " << high_end << " removing " << remove_count << " pixels at each end\n";
result.second = m_feeder.value(high_end);
}
return result;
}
NS_MIA_END
#endif
THistogramFeeder< unsigned char >::size
size_t size() const
THistogramFeeder< unsigned char >::value_type
unsigned char value_type
typedef for generic programming
Definition histogram.hh:104
THistogramFeeder< unsigned char >::THistogramFeeder
THistogramFeeder(unsigned char min, unsigned char max, size_t bins)
Construct the feeder, the parameters are ignored.
THistogramFeeder< unsigned char >::index
size_t index(unsigned char x) const
THistogramFeeder< unsigned char >::value
unsigned char value(size_t k) const
THistogramFeeder
A class to normalize and quantizize input data to a given histogram range with its given number of bi...
Definition histogram.hh:50
THistogramFeeder::value_type
T value_type
typedef for generic programming
Definition histogram.hh:53
THistogramFeeder::size
size_t size() const
Definition histogram.hh:247
THistogramFeeder::THistogramFeeder
THistogramFeeder(T min, T max, size_t bins)
Definition histogram.hh:237
THistogramFeeder::index
size_t index(T x) const
Definition histogram.hh:253
THistogramFeeder::value
T value(size_t k) const
Definition histogram.hh:267
THistogram::median
Feeder::value_type median() const
Definition histogram.hh:368
THistogram::const_iterator
std::vector< size_t >::const_iterator const_iterator
STL iterator.
Definition histogram.hh:139
THistogram::push_range
void push_range(Iterator begin, Iterator end)
Definition histogram.hh:333
THistogram::get_reduced_range
range_type get_reduced_range(double remove) const
Definition histogram.hh:486
THistogram::end
const_iterator end() const
Definition histogram.hh:355
THistogram::THistogram
THistogram(const Feeder &f)
Definition histogram.hh:294
THistogram::value_type
std::pair< typename Feeder::value_type, size_t > value_type
A type for the value-index pair.
Definition histogram.hh:142
THistogram::deviation
double deviation() const
Definition histogram.hh:467
THistogram::average
double average() const
Definition histogram.hh:407
THistogram::push
void push(typename Feeder::value_type x)
Definition histogram.hh:325
THistogram::begin
const_iterator begin() const
Definition histogram.hh:349
THistogram::skewness
double skewness() const
Definition histogram.hh:446
THistogram::excess_kurtosis
double excess_kurtosis() const
Definition histogram.hh:423
THistogram::MAD
Feeder::value_type MAD() const
Definition histogram.hh:381
THistogram::at
const value_type at(size_t idx) const
Definition histogram.hh:398
THistogram::operator[]
size_t operator[](size_t idx) const
Definition histogram.hh:361
THistogram::size
size_t size() const
Definition histogram.hh:319
THistogram::range_type
std::pair< typename Feeder::value_type, typename Feeder::value_type > range_type
Definition histogram.hh:144
defines.hh
NS_MIA_BEGIN
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition defines.hh:33
NS_MIA_END
#define NS_MIA_END
conveniance define to end the mia namespace
Definition defines.hh:36
CUBHistogramFeeder
THistogramFeeder< unsigned char > CUBHistogramFeeder
typedef for the unsigned byte histogram feeder specialization
Definition histogram.hh:120
msgstream.hh
cvdebug
CDebugSink & cvdebug()
Definition msgstream.hh:226