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