histogram.h

/*****************************************************************************
 * VCLib                                                                     *
 * Visual Computing Library                                                  *
 *                                                                           *
 * Copyright(C) 2021-2025                                                    *
 * Visual Computing Lab                                                      *
 * ISTI - Italian National Research Council                                  *
 *                                                                           *
 * All rights reserved.                                                      *
 *                                                                           *
 * This program is free software; you can redistribute it and/or modify      *
 * it under the terms of the Mozilla Public License Version 2.0 as published *
 * by the Mozilla Foundation; either version 2 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              *
 * Mozilla Public License Version 2.0                                        *
 * (https://www.mozilla.org/en-US/MPL/2.0/) for more details.                *
 ****************************************************************************/
 
#ifndef VCL_MATH_HISTOGRAM_H
#define VCL_MATH_HISTOGRAM_H
 
#include <vclib/types.h>
 
#include <algorithm>
#include <cmath>
#include <vector>
 
namespace vcl {
 
// TODO: improve documentation for this class
 
template<typename ScalarType>
class Histogram
{
    std::vector<ScalarType> mHist;   // Counters for bins
    std::vector<ScalarType> mRanges; // Range for bins
 
    ScalarType mMinRange = 0; // Minimum range value
    ScalarType mMaxRange = 1; // Maximum range value
 
    ScalarType mMin = std::numeric_limits<ScalarType>::max();
    ScalarType mMax = std::numeric_limits<ScalarType>::lowest();
 
    // Number of valid bins stored between mMinRange and mMaxRange
    uint mBinNumber = 0;
 
    ScalarType mCnt = 0; // Number of accumulated values
    ScalarType mSum = 0; // Sum of values
    ScalarType mRMS = 0; // Root mean square
 
public:
    Histogram() = default;
 
    Histogram(
        ScalarType minRangeValue,
        ScalarType maxRangeValue,
        uint       nBins,
        ScalarType gamma = 1.0) :
            mMinRange(minRangeValue), mMaxRange(maxRangeValue),
            mBinNumber(nBins)
    {
        mHist.resize(nBins + 2, 0);
 
        mRanges.resize(nBins + 3);
 
        mRanges[0]         = -std::numeric_limits<ScalarType>::max();
        mRanges[nBins + 2] = std::numeric_limits<ScalarType>::max();
 
        double delta = (mMaxRange - mMinRange);
        if (gamma == 1) {
            for (uint i = 0; i <= nBins; ++i)
                mRanges[i + 1] = mMinRange + delta * ScalarType(i) / nBins;
        }
        else {
            for (uint i = 0; i <= nBins; ++i)
                mRanges[i + 1] =
                    mMinRange + delta * std::pow(ScalarType(i) / nBins, gamma);
        }
    }
 
    void clear()
    {
        mHist.clear();
        mRanges.clear();
        mCnt       = 0;
        mSum       = 0;
        mRMS       = 0;
        mBinNumber = 0;
        mMinRange  = 0;
        mMaxRange  = 1;
        mMin       = std::numeric_limits<ScalarType>::max();
        mMax       = std::numeric_limits<ScalarType>::lowest();
    }
 
    void addValue(ScalarType value, ScalarType increment = 1.0)
    {
        uint pos = binIndex(value);
        if (value < mMin)
            mMin = value;
        if (value > mMax)
            mMax = value;
        assert((pos >= 0) && (pos <= mBinNumber + 1));
        mHist[pos] += increment;
        mCnt += increment;
        mSum += value * increment;
        mRMS += (value * value) * increment;
    }
 
    ScalarType minRangeValue() const { return mMinRange; }
 
    ScalarType maxRangeValue() const { return mMaxRange; }
 
    ScalarType sumValues() const { return mSum; }
 
    ScalarType numberValues() const { return mCnt; }
 
    ScalarType minValue() const { return mMin; }
 
    ScalarType maxValue() const { return mMax; }
 
    ScalarType maxBinCount() const
    {
        return *(std::max_element(mHist.begin(), mHist.end()));
    }
 
    ScalarType maxBinCountInRange() const
    {
        return *(std::max_element(mHist.begin() + 1, mHist.end() - 1));
    }
 
    uint binsNumber() const { return mBinNumber; }
 
    ScalarType binCount(uint ind) const { return mHist[ind]; }
 
    ScalarType binLowerBound(uint ind) const { return mRanges[ind]; }
 
    ScalarType binUpperBound(uint ind) const { return mRanges[ind + 1]; }
 
    ScalarType binOfValueCount(ScalarType value) const
    {
        return mHist[binIndex(value)];
    }
 
    ScalarType binOfValueCount(ScalarType value, ScalarType width) const
    {
        return rangeCount(value - width / 2.0, value + width / 2.0);
    }
 
    ScalarType binOfValueWidth(ScalarType value)
    {
        uint pos = BinIndex(value);
        return mRanges[pos + 1] - mRanges[pos];
    }
 
    ScalarType rangeCount(ScalarType rangeMin, ScalarType rangeMax) const
    {
        uint firstBin = binIndex(rangeMin);
        uint lastBin  = binIndex(rangeMax);
 
        ScalarType sum = 0;
        for (uint i = firstBin; i <= lastBin; ++i)
            sum += mHist[i];
        return sum;
    }
 
    ScalarType percentile(ScalarType frac) const
    {
        if (mHist.size() == 0)
            return 0;
 
        // check percentile range
        assert(frac >= 0 && frac <= 1);
 
        ScalarType sum     = mCnt * frac;
        ScalarType partsum = 0;
 
        uint i = 0;
        while (partsum < sum) {
            partsum += mHist[i];
            ++i;
        }
 
        return mRanges[i + 1];
    }
 
    ScalarType average() const { return mSum / mCnt; }
 
    ScalarType rootMeanSquare() const { return std::sqrt(mRMS / mCnt); }
 
    ScalarType variance() const
    {
        return std::abs(mRMS / mCnt - average() * average());
    }
 
    ScalarType standardDeviation() const { return std::sqrt(variance()); }
 
private:
    uint binIndex(ScalarType elem) const
    {
        // 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.
        auto it = std::lower_bound(mRanges.begin(), mRanges.end(), elem);
 
        assert(it != mRanges.begin());
        assert(it != mRanges.end());
 
        uint pos = it - mRanges.begin();
        pos -= 1;
 
        assert(mRanges[pos] < elem);
        assert(elem <= mRanges[pos + 1]);
        return pos;
    }
};
 
/* Specialization Aliases */
 
using Histogramf = Histogram<float>;
using Histogramd = Histogram<double>;
 
} // namespace vcl
 
#endif // VCL_MATH_HISTOGRAM_H