histogram.hxx

/************************************************************************/
/*                                                                      */
/*               Copyright 2011-2012 by Ullrich Koethe                  */
/*                                                                      */
/*    This file is part of the VIGRA computer vision library.           */
/*    The VIGRA Website is                                              */
/*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
/*    Please direct questions, bug reports, and contributions to        */
/*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
/*        vigra@informatik.uni-hamburg.de                               */
/*                                                                      */
/*    Permission is hereby granted, free of charge, to any person       */
/*    obtaining a copy of this software and associated documentation    */
/*    files (the "Software"), to deal in the Software without           */
/*    restriction, including without limitation the rights to use,      */
/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
/*    sell copies of the Software, and to permit persons to whom the    */
/*    Software is furnished to do so, subject to the following          */
/*    conditions:                                                       */
/*                                                                      */
/*    The above copyright notice and this permission notice shall be    */
/*    included in all copies or substantial portions of the             */
/*    Software.                                                         */
/*                                                                      */
/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
/*                                                                      */
/************************************************************************/

#ifndef VIGRA_HISTOGRAM_HXX
#define VIGRA_HISTOGRAM_HXX

#include "config.hxx"
#include "array_vector.hxx"
#include <algorithm>

namespace vigra {

/** \brief Set histogram options.

    HistogramOptions objects are used to pass histogram options to other objects. This \ref acc_hist_options "example" shows how it is is used to pass histogram options to an accumulator chain.
*/
class HistogramOptions
{
  public:

    /** \brief Lower bound for linear range mapping from values to indices. */
    double minimum;

    /** \brief Upper bound for linear range mapping from values to indices. */
    double maximum;

    /** \brief Total number of bins in the histogram. */
    int binCount;

    /** \brief If true, range mapping bounds are defined by minimum and maximum of the data. */
    bool local_auto_init;
    
    /** Initialize members with default values:

    - minimum, maximum = 0.0
    - binCount = 64
    - local_auto_init = false
    */
    HistogramOptions()
    : minimum(0.0), maximum(0.0),
      binCount(64),
      local_auto_init(false)
    {}
    
    /** Set minimum = mi and maximum = ma. Requirement: mi < ma.
    */
    HistogramOptions & setMinMax(double mi, double ma)
    {
        vigra_precondition(mi < ma,
            "HistogramOptions::setMinMax(): min < max required.");
        minimum = mi;
        maximum = ma;
        return *this;
    }

    /** Set binCount = c. Requirement: c > 0.
    */
    HistogramOptions & setBinCount(int c)
    {
        vigra_precondition(c > 0,
            "HistogramOptions::setBinCount(): binCount > 0 required.");
        binCount = c;
        return *this;
    }

    /** Set local_auto_init = true. Requirement: setMinMax() must not have been called before. */
    HistogramOptions & regionAutoInit()
    {
        vigra_precondition(!validMinMax(),
            "HistogramOptions::regionAutoInit(): you must not call setMinMax() when auto initialization is desired.");
        local_auto_init = true;
        return *this;
    }

    /** Set local_auto_init = false. Requirement: setMinMax() must not have been called before. */
    HistogramOptions & globalAutoInit()
    {
        vigra_precondition(!validMinMax(),
            "HistogramOptions::globalAutoInit(): you must not call setMinMax() when auto initialization is desired.");
        local_auto_init = false;
        return *this;
    }
    
    /** Return minimum < maximum.
    */
    bool validMinMax() const
    {
        return minimum < maximum;
    }
};

template <class DataType, class BinType>
class HistogramView
{
    BinType * bins_;
    int size_, stride_;
    DataType offset_;
    double scale_, scaleInverse_;
    
  public:
    HistogramView(DataType const & min, DataType const & max, int binCount, 
                  BinType * bins = 0, int stride = 1)
    : bins_(bins),
      size_(binCount),
      stride_(stride),
      offset_(min),
      scale_(double(binCount) / (max - min)),
      scaleInverse_(1.0 / scale_)
    {}
    
    HistogramView & setData(BinType * bins , int stride = 1)
    {
        bins_ = bins;
        stride_ = stride;
        return *this;
    }
    
    HistogramView & reset()
    {
        if(hasData())
            for(int k=0; k<size_; ++k)
                *(bins_ +k*stride_) = BinType();
        return *this;
    }
    
    void getBinCenters(ArrayVector<DataType> * centers) const
    {
        for(int k=0; k < size_; ++k)
        {
            (*centers)[k] = mapItemInverse(k + 0.5) ;
        }
    }
    
    int size() const
    {
        return size_;
    }
    
    bool hasData() const
    {
        return bins_ != 0;
    }
    
    BinType const & operator[](int k) const
    {
        return *(bins_ + k*stride_);
    }

    double mapItem(DataType const & d) const
    {
        return scale_ * (d - offset_);
    }
    
    DataType mapItemInverse(double d) const
    {
        return DataType(d * scaleInverse_ + offset_);
    }
    
    void add(DataType const & d, BinType weight = NumericTraits<BinType>::one())
    {
        get(int(mapItem(d))) += weight;
    }

  protected:

    BinType & get(int index)
    {
        if(index < 0)
            index = 0;
        if(index >= size_)
            index = size_ - 1;
        return *(bins_ + index*stride_);
    }
};

template <class T>
class TrapezoidKernel
{
  public:
    typedef T value_type;
    
    T operator[](double f) const
    {
        if(f < -0.5)
            return 0.5*(f + 1.5);
        if(f > 0.5)
            return 0.5*(1.5 - f);
        return 0.5;
    }
    
    double radius() const
    {
        return 1.5;
    }
    
    T findMaximum(double l, double c, double r) const
    {
        double curv = -2.0*c + r + l;
        if(curv == 0.0)
            return T(-0.5);
        double extr = 0.5*(l-r) / curv;
        if(curv < 0.0)
        {
            return extr < -0.5
                       ? T(-0.5)
                       : extr > 0.5
                              ? T(0.5)
                              : T(extr);
        }
        else
        {
            return extr < 0.0
                       ? T(0.5)
                       : T(-0.5);
        }
    }
    
    bool findMode(double l, double c, double r, double * m) const
    {
        double curv = -2.0*c + r + l;
        if(curv >= 0.0)
            return false;
        *m = 0.5*(l-r) / curv;
        if(*m < -0.5 || *m > 0.5)
            return false;
        return true;
    }
};

template <class DataType, class KernelType>
class KernelHistogramView
: public HistogramView<DataType, typename KernelType::value_type>
{
    KernelType kernel_;
    int radius_;
    
  public:
  
    typedef typename KernelType::value_type BinType;
    typedef HistogramView<DataType, BinType> BaseType;

    KernelHistogramView(DataType const & min, DataType const & max, int binCount, 
                        BinType * bins = 0, int stride = 1)
    : BaseType(min, max, binCount, bins, stride),
      radius_(kernel_.radius()-0.5) // FIXME: this needs generalization
    {}
    
    void add(DataType const & d, BinType weight = NumericTraits<BinType>::one())
    {
        double mapped = this->mapItem(d);
        double f = mapped - std::floor(mapped) - kernel_.radius();
        int center = int(mapped);
        
        for(int k=center+radius_; k>=center-radius_; --k, f += 1.0)
        {   
            this->get(k) += weight*kernel_[f];
        }
    }
    
    DataType findMode() const
    {
        double mmax = 0, vmax = 0, m;
        
        for(int k=0; k<this->size(); ++k)
        {
            double l = k > 0
                         ? (*this)[k-1]
                         : 0.0;
            double c = (*this)[k];
            double r = k < this->size() - 1
                         ? (*this)[k+1]
                         : 0.0;
            if(kernel_.findMode(l, c, r, &m))
            {
                double v = l*kernel_[m+1.0] + c*kernel_[m] + r*kernel_[m-1.0];
                if(vmax < v)
                {
                    mmax = m + k + 0.5;
                    vmax = v;
                }
            }
        }
        return this->mapItemInverse(mmax);
    }
    
    template <class Array>
    void findModes(Array * modes)
    {
        double m;
        for(int k=0; k<this->size(); ++k)
        {
            double l = k > 0
                         ? (*this)[k-1]
                         : 0.0;
            double c = (*this)[k];
            double r = k < this->size() - 1
                         ? (*this)[k+1]
                         : 0.0;
            if(kernel_.findMode(l, c, r, &m))
            {
                double v = l*kernel_[m+1.0] + c*kernel_[m] + r*kernel_[m-1.0];
                modes->push_back(std::make_pair(this->mapItemInverse(m + k + 0.5), v));
            }
        }
    }
};

template <class DataType, class BinType>
class Histogram
: public HistogramView<DataType, BinType>
{
  public:
    typedef HistogramView<DataType, BinType> BaseType;
    ArrayVector<BinType> data_;
    
  public:
    Histogram(DataType const & min, DataType const & max, int binCount, 
                  BinType * /*bins*/ = 0, int /*stride*/ = 1)
    : BaseType(min, max, binCount),
      data_(binCount)
    {
        this->setData(&data_[0]);
    }
    
    Histogram const & reset()
    {
        this->setData(&data_[0]);
        BaseType::reset();
        return *this;
    }
 };

} // namespace vigra

#endif // VIGRA_HISTOGRAM_HXX