CRNHistogram.cpp

Go to the documentation of this file.

/* Copyright 2006-2016 Yann LEYDIER, INSA-Lyon, CoReNum, Université Paris Descartes, ENS-Lyon
 *
 * This file is part of libcrn.
 *
 * libcrn is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * libcrn 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 Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with libcrn.  If not, see <http://www.gnu.org/licenses/>.
 *
 * file: SHistogram cpp
 * \author Yann LEYDIER, Jean DUONG, Asma OUJI, Frank LeBOURGEOIS
 */

#include <CRNStatistics/CRNHistogram.h>
#include <CRNData/CRNData.h>
#include <CRNException.h>
#include <CRNImage/CRNImageBW.h>
#include <CRNData/CRNDataFactory.h>
#include <CRNIO/CRNIO.h>
#include <functional>
#include <CRNi18n.h>

using namespace crn;

/*****************************************************************************/
Histogram::Histogram(size_t s, unsigned int v):
    bins(s, v),
    compression(1)
{ }

/*****************************************************************************/
Histogram::Histogram(const Histogram &src, unsigned int c):
    compression(c)
{

    if (c == 1)
    {
        bins = src.bins;
    }
    else
    {
        size_t q = src.Size() / compression;
        size_t r = src.Size() % compression;
        size_t newsize;
        if (r > 0)
        {
            newsize = q + 1;
        }
        else
        {
            newsize = q;
        }

        datatype(newsize, 0).swap(bins);
        size_t k = 0;

        for (size_t i = 0; i < q; i++)
        {
            unsigned int sum = 0;

            for (unsigned int j = 0; j < c; j++)
            {
                sum += src.bins[k];
                k++;
            }

            bins[i] = sum;
        }

        if (r > 0)
        {
            unsigned int sum = 0;

            for (size_t j = k; j < src.Size(); j++)
            {
                sum += src.bins[k];
                k++;
            }

            bins[q] = sum;
        }
    }
}

/*****************************************************************************/
void Histogram::SetBin(size_t k, unsigned int v)
{
    if (k >= bins.size())
    {
        throw ExceptionDomain(StringUTF8("Histogram::SetBin(size_t k, unsigned int v): ") +
                _("invalid index."));
    }

    bins[k] = v;
}

/*****************************************************************************/
void Histogram::IncBin(size_t k, unsigned int i)
{
    if (k >= bins.size())
    {
        throw ExceptionDomain(StringUTF8("Histogram::IncBin(size_t k, unsigned int i): ") +
                _("invalid index."));
    }

    bins[k] += i;
}

/*****************************************************************************/
unsigned int Histogram::GetBin(size_t k) const
{
    if (k >= bins.size())
    {
        throw ExceptionDomain(StringUTF8("Histogram::GetBin(size_t k): ") +
                _("invalid index."));
    }

    return bins[k];
}


/*****************************************************************************/
unsigned int Histogram::CumulateBins() const
{
    unsigned int Cumul = 0;

    for (auto & elem : bins)
    {
        Cumul += elem;
    }

    return Cumul;
}

/*****************************************************************************/
double Histogram::Mean() const
{
    double m = 0.0;
    unsigned int p = 0;

    for (size_t k = 0; k < bins.size(); k++)
    {
        m += double(k * bins[k]);
        p += bins[k];
    }

    m /= p;

    return m;
}

/*****************************************************************************/
double Histogram::Variance(double m) const
{
    double v = 0.0;
    unsigned int p=0;

    for (size_t k = 0; k < bins.size(); k++)
    {
        double d = double(k) - m;

        v += d * d * bins[k];
        p += bins[k];
    }

    v /= p;

    return v;
}

/*****************************************************************************/
double Histogram::Variance() const
{
    double m = Mean();
    double v = 0.0;
    unsigned int p=0;

    for (size_t k = 0; k < bins.size(); k++)
    {
        double d = double(k) - m;

        v += d * d * bins[k];
        p += bins[k];
    }

    v /= p;

    return v;
}

/*****************************************************************************/
unsigned int Histogram::Max() const
{
    unsigned int m = 0;

    for (auto & elem : bins)
    {
        m = crn::Max(m, elem);
    }

    return m;
}

/*****************************************************************************/
unsigned int Histogram::Min() const
{
    unsigned int m = bins[0];

    for (auto & elem : bins)
    {
        m = crn::Min(m, elem);
    }

    return m;
}

/*****************************************************************************/
size_t Histogram::Argmax() const
{
    unsigned int m = 0;
    size_t a = 0;

    for (size_t k = 0; k < bins.size(); k++)
    {
        unsigned int v = bins[k];

        if (v > m)
        {
            m = bins[k];
            a = k;
        }
    }

    return a;
}

/*****************************************************************************/
size_t Histogram::Argmin() const
{
    unsigned int m = bins[0];
    size_t a = 0;

    for (size_t k = 0; k < bins.size(); k++)
    {
        unsigned int v = bins[k];

        if (v < m)
        {
            m = bins[k];
            a = k;
        }
    }

    return a;
}

/*****************************************************************************/
void Histogram::SetCeiling(unsigned int m)
{
    for (auto & elem : bins)
    {
        elem = crn::Min(m,elem);

    }
}

/*****************************************************************************/
void Histogram::ScaleMaxTo(unsigned int m)
{
    unsigned int c = Max(); /* Maximal value for this histogram */

    if (c == 0)
    {
        CRNWarning(String(U"Histogram::ScaleMaxTo(int m): ") +
                _("flat histogram."));
        return;
    }

    double s = (double) m / (double) c; /* Scale factor */

    for (auto & elem : bins)
    {
        unsigned int v = (unsigned int)(elem * s);

        elem = v;
    }
}

/*****************************************************************************/
void Histogram::AverageSmoothing(size_t d)
{
    datatype H(bins.size());

    for (size_t k = 0; k < bins.size(); k++)
    {
        int left = crn::Max(0, int(k) - int(d));
        int right = crn::Min(int(bins.size()) - 1, int(k) + int(d));

        uint64_t m = 0;
        uint64_t n = right - left + 1;

        for (int j = left; j <= right; j++)
        {
            m += bins[j];
        }

        H[k] = (unsigned int)(m / n);
    }

    bins.swap(H);
}

/*****************************************************************************/
void Histogram::CircularAverageSmoothing(size_t d)
{
    if (d > bins.size())
    {
        CRNWarning(String(U"Histogram::CircularAverageSmoothing(size_t d): ") +
                _("window larger than histogram. Cropping."));
        d = bins.size(); // wew, dirty !
    }

    datatype H(bins.size());

    for (size_t k = 0; k < bins.size(); k++)
    {
        int left = int(k) - int(d);
        int right = int(k + d);
        double m = 0.0;
        double n = double(right) - double(left) + 1;

        for (int j = left; j <= right; j++)
        {
            int index = j;
            if (index < 0) index = int(bins.size()) + index;
            if (index >= int(bins.size())) index -= int(bins.size());
            m += (double) bins[index] / n;
        }

        H[k]=(unsigned int)m;
    }

    bins.swap(H);
}

/*****************************************************************************/
std::vector<size_t> Histogram::Modes() const
{
    std::vector<size_t> modes;
    if (bins.empty())
        return modes;
    if (bins.size() == 1)
    {
        modes.push_back(0);
        return modes;
    }
    // first value is a mode ?
    if (bins[0] > bins[1])
        modes.push_back(0);
    for (size_t b = 1; b < bins.size() - 1; ++b)
    {
        if (bins[b] > bins[b - 1])
        { // increasing
            size_t first = b;
            if (bins[b] < bins[b + 1])
            { // still increasing, pass
                continue;
            }
            while (bins[b] == bins[b + 1])
            { // plateau
                b += 1;
                if (b == bins.size() - 1)
                { // reached the end of the histogram
                    modes.push_back((first + bins.size() - 1) / 2);
                    break;
                }
                else if (b == bins.size() - 2)
                { // reached the penultimate
                    if (bins[b] == bins[b + 1])
                    { // plateau till the end
                        modes.push_back((first + bins.size() - 1) / 2);
                    }
                    if (bins[b] > bins[b + 1])
                    { // peak
                        modes.push_back((first + b) / 2);
                    }
                    break;
                }
            }
            if ((b < bins.size() - 1) && (bins[b] > bins[b + 1]))
            { // peak !
                modes.push_back((first + b) / 2);
            }
        }
    }
    if (bins[bins.size() - 2] < bins[bins.size() - 1])
        modes.push_back(bins.size() - 1);

    return modes;
}

std::vector<size_t> Histogram::StableModes() const
{
    // copy histogram
    Histogram s(*this);

    std::map<size_t, size_t> mcount; // key = number of modes, value = number of occurrences
    std::map<size_t, std::vector<size_t> > msum; // key = number of modes, value = sum of the modes

    std::vector<size_t> modes(s.Modes());
    size_t nmodes = modes.size();
    size_t niter = 0;
    while (nmodes > 1)
    {
        auto cit = mcount.find(nmodes);
        if (cit != mcount.end())
        { // update
            cit->second += 1;
            auto sit = msum.find(nmodes);
            std::transform(sit->second.begin(), sit->second.end(), modes.begin(), sit->second.begin(), std::plus<size_t>());
        }
        else
        { // add
            mcount.insert(std::make_pair(nmodes, 1));
            msum.insert(std::make_pair(nmodes, modes));
        }
        // smooth
        s.AverageSmoothing(1);
        modes = s.Modes();
        nmodes = modes.size();
        niter += 1;
    }

    // if there is only one iteration, then it means there is only one mode
    if (niter == 1)
    {
        return modes;
    }

    // look for most stable mode number
    nmodes = 0;
    size_t maxpop = 0;
    for (auto & elem : mcount)
    {
        if (elem.second > maxpop)
        {
            maxpop = elem.second;
            nmodes = elem.first;
        }
    }

    // compute the vector of modes
    if (nmodes == 0)
        return std::vector<size_t>();

    modes.swap(msum[nmodes]);
    std::transform(modes.begin(), modes.end(), modes.begin(), std::bind2nd(std::divides<size_t>(), maxpop));
    return modes;
}

/*****************************************************************************/
Histogram Histogram::MakeIntersection(const Histogram &h) const
{
    if (bins.size() != h.Size())
    {
        throw ExceptionDimension(StringUTF8("Histogram::Intersection(const Histogram &h): ") +
                _("histograms must have same size."));
    }

    Histogram i(bins.size());

    for (size_t k = 0; k < bins.size(); k++)
    {
        i.bins[k] = (unsigned int)(crn::Min(h.GetBin(k),bins[k]));
    }

    return i;
}

/*****************************************************************************/
double Histogram::IntersectionDivergence(const Histogram &h) const
{
    if (bins.size() != h.Size())
    {
        throw ExceptionDomain(StringUTF8("Histogram::IntersectionDivergence("
                    "const Histogram &h): ") + _("histograms must have same size."));
    }

    double cumulMin = 0.0;
    double cumulH = h.CumulateBins();

    for (size_t k = 0; k < bins.size(); k++)
    {
        cumulMin += crn::Min(h.GetBin(k),bins[k]);
    }

    return (1.0 - (cumulMin / cumulH));
}

/*****************************************************************************/
double Histogram::Correlation(const Histogram &h) const
{
    if (bins.size() != h.Size())
    {
        throw ExceptionDimension(StringUTF8("Histogram::Correlation("
                        "const Histogram &h): ") + _("histograms must have same size."));
    }

    int offset_1 = (int)(CumulateBins() / bins.size());
    int offset_2 = (int)(h.CumulateBins() / bins.size());

    double cumul_numer = 0.0;
    double cumul_denom_1 = 0.0;
    double cumul_denom_2 = 0.0;

    for (size_t k = 0; k < bins.size(); k++)
    {
        int v_1 = int(bins[k]) - offset_1;
        int v_2 = int(h.bins[k]) - offset_2;

        cumul_numer += v_1 * v_2;
        cumul_denom_1 += v_1 * v_1;
        cumul_denom_2 += v_2 * v_2;
    }

    cumul_numer /= sqrt(cumul_denom_1);
    cumul_numer /= sqrt(cumul_denom_2);

    return cumul_numer;
}


/*****************************************************************************/
double Histogram::Chi2(const Histogram &h) const
{
    if (bins.size() != h.Size())
    {
        throw ExceptionDimension(StringUTF8("Histogram::IntersectionDivergence("
                        "const Histogram &h): ") + _("histograms must have same size."));
        return HUGE_VAL;
    }

    double cumul_numer = 0.0;
    double cumul_denom = 0.0;

    for (size_t k = 0; k < bins.size(); k++)
    {
        auto v_1 = int(bins[k]);
        auto v_2 = int(h.bins[k]);
        int diff = v_1 - v_2;

        cumul_numer += diff * diff;
        cumul_denom += v_1 + v_2;
    }

    return cumul_numer / cumul_denom;
}

/*****************************************************************************/
double Histogram::MinkowskiDistance(const Histogram &h, double r) const
{
    if (bins.size() != h.Size())
    {
        throw ExceptionDimension(StringUTF8("Histogram::MinkowskiDistance("
                    "const Histogram &h, double r): ") +
                _("histograms must have same size."));
    }

    double cumul = 0.0;
    double invR = 1.0 / (double)(r);

    for (size_t k = 0; k < bins.size(); k++)
    {
        cumul += pow(fabs((double)(h.GetBin(k)) - (double)(bins[k])),r);
    }

    return pow(cumul, invR);
}

/*****************************************************************************/
double Histogram::JeffreyDivergence(const Histogram &h) const
{
    if (bins.size() != h.Size())
    {
        throw ExceptionDimension(StringUTF8("Histogram::JeffreyDistance(const Histogram &h): ") +
                _("histograms must have same size."));
    }

    double cumul = 0.0;

    for (size_t k = 0; k < bins.size(); k++)
    {
        double h1 = bins[k];
        double h2 = h.GetBin(k);
        double m = (h1 + h2) / 2.0;

        cumul += h1 * log(h1 / m) + h2 * log(h2 / m);
    }

    return cumul;
}

/*****************************************************************************/
double Histogram::MatchDistance(const Histogram &h) const
{
    if (bins.size() != h.Size())
    {
        throw ExceptionDimension(StringUTF8("Histogram::MatchDistance(const Histogram &h): ") +
                _("histograms must have same size."));
    }

    double cumul = 0.0;
    double cumulH1  = 0.0;
    double cumulH2  = 0.0;

    for (size_t k = 0; k < bins.size(); k++)
    {
        cumulH1 += bins[k];
        cumulH2 += h.GetBin(k);
        cumul += fabs(cumulH1 - cumulH2);
    }

    return cumul;
}

/*****************************************************************************/
double Histogram::KolmogorovSmirnovDistance(const Histogram &h) const
{
    if (bins.size() != h.Size())
    {
        throw ExceptionDimension(StringUTF8("Histogram::KolmogorovSmirnovDistance("
                    "const Histogram &h): ")+ _("histograms must have same size."));
    }

    double max = 0.0;
    double cumulH1  = 0.0;
    double cumulH2  = 0.0;

    for (size_t k = 0; k < bins.size(); k++)
    {
        cumulH1 += bins[k];
        cumulH2 += h.GetBin(k);
        double diff = fabs(cumulH1 - cumulH2);

        if (diff > max)
        {
            max = diff;
        }
    }

    return max;
}

double Histogram::EMD(const Histogram &h) const
{
    if (bins.size() != h.Size())
    {
        throw ExceptionDimension(StringUTF8("Histogram::EMD(const Histogram &h): ")+ _("histograms must have same size."));
    }
    auto h1 = *this;
    h1.Cumulate();
    auto h2 = h;
    h2.Cumulate();
    auto dist = 0.0;
    for (auto tmp : crn::Range(h1))
        dist += double(Abs(int(h1[tmp]) - int(h2[tmp])));
    return dist / double(h1.Size());
}

double Histogram::CEMD(const Histogram &h) const
{
    if (bins.size() != h.Size())
    {
        throw ExceptionDimension(StringUTF8("Histogram::CEMD(const Histogram &h): ")+ _("histograms must have same size."));
    }
    auto dist = std::numeric_limits<double>::max();
    // compute EMD for each possible cumulative histogram
    for (auto k : crn::Range(h))
    {
        // compute cumulative histogram
        auto h1 = Histogram(Size());
        auto h2 = Histogram(Size());
        // sum from k to end
        auto sum1 = 0u;
        auto sum2 = 0u;
        for (auto i = k; i < Size(); ++i)
        {
            sum1 += (*this)[i];
            sum2 += h[i];
        }
        // cumulate from k to end and 0 to i (i < k)
        for (auto i = size_t(0); i < k; ++i)
        {
            h1[i] = sum1;
            h2[i] = sum2;
            for (auto j = size_t(0); j <= i; ++j)
            {
                h1[i] += (*this)[j];
                h2[i] += h[j];
            }
        }
        // cumulate from k to i (i >= k)
        for (auto i = k; i < Size(); ++i)
            for (auto j = k; j <= i; ++j)
            {
                h1[i] += (*this)[j];
                h2[i] += h[j];
            }
        // compute EMPD
        auto d = 0.0;
        for (auto tmp : crn::Range(h1))
            d += double(Abs(int(h1[tmp]) - int(h2[tmp])));
        // keep min
        if (d < dist)
            dist = d;
    }
    return dist / double(Size());
}

/*****************************************************************************/
Histogram::~Histogram()
{ }

/*****************************************************************************/
void Histogram::Append(const Histogram &h)
{
    bins.insert(bins.end(), h.bins.begin(), h.bins.end());
}

/*****************************************************************************/
void Histogram::Resize(size_t newsize)
{
    if (newsize == 0)
        throw crn::ExceptionDomain(_("Cannot resize histogram with null size."));

    if (newsize == bins.size())
        return;

    datatype newbins(newsize);

    double ratio = (double)bins.size() / (double)newsize;
    double begin = 0;
    for (size_t tmp = 0; tmp < newsize; tmp++)
    {
        double val = 0;
        auto end = double(tmp + 1) * ratio;
        if (end >= double(bins.size())) end = double(bins.size()) - 0.0001; // yay... -_-
        double tb = begin;
        while (floor(tb) != floor(end))
        {
            double te = ceil(tb);
            if (te == tb) te += 1;
            val += (te - tb) * bins[(size_t)floor(tb)];
            tb = te;
        }
        val += (end - tb) * bins[(size_t)floor(tb)];
        if ((end - begin) != 0)
            val /= (end - begin);
        //GET: using round in place of direct cast
        //GETnewbins[tmp] = (unsigned int)val;
        newbins[tmp] = (unsigned int)round(val);
        begin = end;
    }

    bins.swap(newbins);
}

/*****************************************************************************/
void Histogram::Deserialize(xml::Element &el)
{
    bins.clear();
    if (el.GetValue() != "Histogram")
    {
        throw ExceptionInvalidArgument(StringUTF8("bool Histogram::deserialize(xml::Element &el): ") +
                _("Wrong XML element."));
    }
    xml::Node n(el.GetFirstChild());
    if (!n)
    {
        throw ExceptionNotFound(StringUTF8("bool Histogram::deserialize(xml::Element &el): ") +
                _("Cannot get CDATA."));
    }
    xml::Text t = n.AsText(); // may throw
    auto v = Data::ASCII85Decode<unsigned int>(t.GetValue());
    if (v.empty())
    {
        throw ExceptionRuntime(StringUTF8("bool Histogram::deserialize(xml::Element &el): ") +
                _("Cannot convert CDATA."));
    }
    bins.swap(v);
}

/*****************************************************************************/
xml::Element Histogram::Serialize(xml::Element &parent) const
{
    xml::Element el(parent.PushBackElement("Histogram"));
    auto text = Data::ASCII85Encode(reinterpret_cast<const uint8_t * const>(bins.data()), bins.size() * sizeof(datatype::value_type));
    xml::Text t(el.PushBackText(text, false));

    return el;
}

String Histogram::ToString() const
{
    String s(U"Histogram: ");
    for (auto & elem : bins)
    {
        s += elem;
        s += U" ";
    }
    return s;
}

/*****************************************************************************/
ImageBW Histogram::MakeImageBW(size_t height) const
{
    auto bw = ImageBW(bins.size(), height, pixel::BWWhite);
    auto val = double(height) / double(Max());
    for (size_t k = 0; k < bins.size(); k++)
    {
        auto bin = GetBin(k);
        for (int i = 0; i < int(double(bin) * val) - 1; ++i)
        {
            bw.At(k, height - 1 - i) = pixel::BWBlack;
        }
    }
    return bw;
}

void Histogram::Cumulate()
{
    unsigned int sum = 0;
    for (size_t k = 0; k < bins.size(); k++)
    {
        sum += GetBin(k);
        SetBin(k, sum);
    }
}

/*****************************************************************************/
ImageBW Histogram::MakeRadialImageBW(size_t radius) const
{
    Rect r;
    double factor = crn::Twice(M_PI) / double(bins.size());
    double max = Max();
    for (size_t tmp = 0; tmp < bins.size(); tmp++)
    {
        double tx = double(GetBin(tmp) * radius) / max * cos(double(tmp) * factor);
        double ty = double(GetBin(tmp) * radius) / max * sin(double(tmp) * factor);
        r |= Rect((int)tx, (int)ty);
    }
    if (!r.GetWidth() || !r.GetHeight())
        return ImageBW{};
    ImageBW img = ImageBW(r.GetWidth(), r.GetHeight(), pixel::BWWhite);
    for (size_t tmp = 0; tmp < bins.size(); tmp++)
    {
        double tx = double(GetBin(tmp) * radius) / max * cos(double(tmp) * factor);
        double ty = double(GetBin(tmp) * radius) / max * sin(double(tmp) * factor);
        img.DrawLine(-r.GetLeft(), -r.GetTop(), (int)tx - r.GetLeft(),
                (int)ty - r.GetTop(), pixel::BWBlack);
    }
    return img;
}

/*****************************************************************************/
size_t Histogram::Fisher() const
{
    std::vector<bool> notempty(bins.size());
    std::vector<double> tab(bins.size());

    for (size_t tmp = 0; tmp < bins.size(); tmp++)
    {
        if (bins[tmp] > 0)
        {
            notempty[tmp] = true;
            tab[tmp] = double((tmp + 1) * bins[tmp]);
        }
        else
        {
            notempty[tmp] = false;
            tab[tmp] = 0;
        }
    }

    double s1 = 0, s2 = 0;
    long n1 = 0, n2 = 0;
    for (size_t tmp = 0; tmp < bins.size(); tmp++)
    {
        if (notempty[tmp])
        {
            s2 += tab[tmp];
            n2 += bins[tmp];
        }
    }
    size_t index = 0;
    double f(0.0);
    do
    {
        if (notempty[index])
        {
            s1 += tab[index];
            n1 += bins[index];
            s2 -= tab[index];
            n2 -= bins[index];
        }
        if ((n1 * n2) == 0)
            f = -20.0;
        else
            f = double(index) - ((s1 / double(n1)) + (s2 / double(n2))) / 2.0;
        index += 1;
    } while ((index < bins.size()) && (f <= 0.0));

    if (index)
        return index - 1;
    else
        return 0;
}

/*****************************************************************************/
size_t Histogram::EntropyThreshold() const
{
    std::vector<double> entropy(bins.size());
    for (size_t tmp = 0; tmp < bins.size(); tmp++)
    {
        if (bins[tmp])
            entropy[tmp] = bins[tmp] * log((double)bins[tmp]);
        else
            entropy[tmp] = 0;
    }
    double max = 0;
    size_t thresh = 0;
    for (size_t t = 0; t < bins.size(); t++)
    {
        double s1 = 0, s2 = 0, n1 = 0, n2 = 0;
        for (size_t tmp = 0; tmp <= t; tmp++)
            if (bins[tmp])
            {
                s1 += entropy[tmp];
                n1 += bins[tmp];
            }
        for (size_t tmp = t + 1; tmp < bins.size(); tmp++)
            if (bins[tmp])
            {
                s2 += entropy[tmp];
                n2 += bins[tmp];
            }
        double en = 0;
        if ((n1 * n2) != 0)
            en = -(s1/n1) - (s2/n2) + log(n1*n2);
        if (en > max)
        {
            max = en;
            thresh = t;
        }
    }
    return thresh;
}

/*****************************************************************************/

double Histogram::Entropy() const
{
    auto e = 0.0;
    auto c = double(CumulateBins());

    for (auto & elem : bins)
    {
        auto b_k = double(elem);

        if (b_k != 0)
        {
            auto p_k = double(b_k / c);

            e += p_k * log(p_k);
        }
    }

    return -e;
}

/*****************************************************************************/

size_t Histogram::MedianValue() const
{
    if (bins.empty())
        throw ExceptionUninitialized(_("The histogram is empty."));
    unsigned int mid = CumulateBins() / 2;
    unsigned int sum = 0;
    size_t index;
    for (index = 0; index < bins.size(); index++)
    {
        sum += bins[index];
        if (sum >= mid)
            break;
    }
    return index;
}

Histogram Histogram::MakePopulationHistogram() const
{
    Histogram poph(Max() + 1);
    for (auto & elem : bins)
    {
        poph.IncBin(elem);
    }
    return poph;
}

CRN_BEGIN_CLASS_CONSTRUCTOR(Histogram)
    CRN_DATA_FACTORY_REGISTER(U"Histogram", Histogram)
    Cloner::Register<Histogram>();
    Ruler::Register<Histogram>();
CRN_END_CLASS_CONSTRUCTOR(Histogram)