多级 (4) Otsu 阈值

Question

我正在尝试实现多级 Otsu 的阈值，更具体地说，我需要 3 个阈值/4 个类。

我知道关于它的 2 个类似问题：#34856019 和 #22706742。 问题是我没有得到好的结果：我已经阅读了几篇文章，其中包含示例图像，并且该代码发现的阈值与这些论文中的不同。

假设我有一张黑色背景上有 3 个圆圈的图片，圆圈的亮度从非常亮到暗不等：

Sample Image

我是否正确假设得到 4 个类别：黑色背景和 3 个根据圆圈强度的类别？

我的程序给了我这些阈值：226、178、68

因此，第三个圆圈完全不可见 - 它与背景属于同一类。

有人可以检查这些值和/或源代码吗？也许可以使用 Matlab 或其他方式检查此图像... 顺便说一句，处理除以零的最佳方法是什么，这通常发生在直方图中的零值？ 源代码：

void MultilevelThresholding(cv::Mat& src)
{
    int histogram[256] = { 0 };
    int pixelsCount = src.cols * src.rows;

    for (int y = 0; y < src.rows; y++)
    {
        for (int x = 0; x < src.cols; x++)
        {
            uchar value = src.at<uchar>(y, x);
            histogram[value]++;
        }
    }

    double c = 0;
    double Mt = 0;

    double p[256] = { 0 };
    for (int i = 0; i < 256; i++)
    {
        p[i] = (double) histogram[i] / (double) pixelsCount;
        Mt += i * p[i];
    }

    int optimalTreshold1 = 0;
    int optimalTreshold2 = 0;
    int optimalTreshold3 = 0;

    double maxBetweenVar = 0;

    double w0 = 0;
    double m0 = 0;
    double c0 = 0;
    double p0 = 0;

    double w1 = 0;
    double m1 = 0;
    double c1 = 0;
    double p1 = 0;

    double w2 = 0;
    double m2 = 0;
    double c2 = 0;
    double p2 = 0;
    for (int tr1 = 0; tr1 < 256; tr1++)
    {
        p0 += p[tr1];
        w0 += (tr1 * p[tr1]);
        if (p0 != 0)
        {
            m0 = w0 / p0;
        }

        c0 = p0 * (m0 - Mt) * (m0 - Mt);

        c1 = 0;
        w1 = 0;
        m1 = 0;
        p1 = 0;
        for (int tr2 = tr1 + 1; tr2 < 256; tr2++)
        {

            p1 += p[tr2];
            w1 += (tr2 * p[tr2]);
            if (p1 != 0)
            {
                m1 = w1 / p1;
            }

            c1 = p1 * (m1 - Mt) * (m1 - Mt);


            c2 = 0;
            w2 = 0;
            m2 = 0;
            p2 = 0;
            for (int tr3 = tr2 + 1; tr3 < 256; tr3++)
            {

                p2 += p[tr3];
                w2 += (tr3 * p[tr3]);
                if (p2 != 0)
                {
                    m2 = w2 / p2;
                }

                c2 = p2 * (m2 - Mt) * (m2 - Mt);

                c = c0 + c1 + c2;

                if (maxBetweenVar < c)
                {
                    maxBetweenVar = c;
                    optimalTreshold1 = tr1;
                    optimalTreshold2 = tr2;
                    optimalTreshold3 = tr3;
                }
            }
        }
    }

Answer 1

所以，我想通了。 4类（3个阈值）Otsu阈值的最终源代码：

// cv::Mat& src - source image's matrix
    int histogram[256] = { 0 };
    int pixelsCount = src.cols * src.rows;

    for (int y = 0; y < src.rows; y++)
    {
        for (int x = 0; x < src.cols; x++)
        {
            uchar value = src.at<uchar>(y, x);
            histogram[value]++;
        }
    }

    double c = 0;
    double Mt = 0;

    double p[256] = { 0 };
    for (int i = 0; i < 256; i++)
    {
        p[i] = (double) histogram[i] / (double) pixelsCount;
        Mt += i * p[i];
    }

    int optimalTreshold1 = 0;
    int optimalTreshold2 = 0;
    int optimalTreshold3 = 0;

    double maxBetweenVar = 0;

    double w0 = 0;
    double m0 = 0;
    double c0 = 0;
    double p0 = 0;

    double w1 = 0;
    double m1 = 0;
    double c1 = 0;
    double p1 = 0;

    double w2 = 0;
    double m2 = 0;
    double c2 = 0;
    double p2 = 0;
    for (int tr1 = 0; tr1 < 256; tr1++)
    {
        p0 += p[tr1];
        w0 += (tr1 * p[tr1]);
        if (p0 != 0)
        {
            m0 = w0 / p0;
        }

        c0 = p0 * (m0 - Mt) * (m0 - Mt);

        c1 = 0;
        w1 = 0;
        m1 = 0;
        p1 = 0;
        for (int tr2 = tr1 + 1; tr2 < 256; tr2++)
        {

            p1 += p[tr2];
            w1 += (tr2 * p[tr2]);
            if (p1 != 0)
            {
                m1 = w1 / p1;
            }

            c1 = p1 * (m1 - Mt) * (m1 - Mt);


            c2 = 0;
            w2 = 0;
            m2 = 0;
            p2 = 0;
            for (int tr3 = tr2 + 1; tr3 < 256; tr3++)
            {

                p2 += p[tr3];
                w2 += (tr3 * p[tr3]);
                if (p2 != 0)
                {
                    m2 = w2 / p2;
                }

                c2 = p2 * (m2 - Mt) * (m2 - Mt);

                double p3 = 1 - (p0 + p1 + p2);
                double w3 = Mt - (w0 + w1 + w2);
                double m3 = w3 / p3;
                double c3 = p3 * (m3 - Mt) * (m3 - Mt);

                double c = c0 + c1 + c2 + c3;

                if (maxBetweenVar < c)
                {
                    maxBetweenVar = c;
                    optimalTreshold1 = tr1;
                    optimalTreshold2 = tr2;
                    optimalTreshold3 = tr3;
                }
            }
        }
    }

Source image

Result: 3 thresholds / 4 classes

阈值：179、92、25