多级阈值。大津的方法答案

【问题标题】：Multi-Level Thresholding. Otsu's Method多级阈值。大津的方法
【发布时间】：2016-08-30 06:31:54
【问题描述】：

我知道这个问题可能与Multi otsu(multi-thresholding) with openCV 重复，因为它非常相似。但我没有足够的权限对此发表评论，所以我的问题得到了回答。

我正在尝试使用 Otsu 使用 OPENCV 库的方法来做三级阈值。我按照@Antoni4 发布的指南进行操作，但它仅涵盖两级阈值，我无法将其扩展到三级。

我知道原始 Otsu 的方法如何工作，最大化类间方差，但我不知道如何添加外部（或内部）for 循环来计算第三级阈值。

我正在将此代码用于两级阈值：

#include "Shadow01-1.cuh"

void multiThresh(double &optimalThresh1, double &optimalThresh2,double &optimalThresh3, cv::Mat &imgHist, cv::Mat &src) {
double W0K, W1K, W2K, W3K, M0, M1, M2, M3, currVarB, maxBetweenVar, M0K, M1K, M2K, M3K, MT;

unsigned char *histogram = (unsigned char*)(imgHist.data);

int N = src.rows*src.cols;
W0K = 0;
W1K = 0;
M0K = 0;
M1K = 0;
MT = 0;
maxBetweenVar = 0;

for (int k = 0; k <= 255; k++) {
    MT += k * (histogram[k] / (double) N);
}

for (int t1 = 0; t1 <= 255; t1++)
{
    W0K += histogram[t1] / (double) N; //Pi
    M0K += t1 * (histogram[t1] / (double) N); //i * Pi
    M0 = M0K / W0K; //(i * Pi)/Pi

    W1K = 0;
    M1K = 0;

    for (int t2 = t1 + 1; t2 <= 255; t2++)
    {
        W1K += histogram[t2] / (double) N; //Pi
        M1K += t2 * (histogram[t2] / (double) N); //i * Pi
        M1 = M1K / W1K; //(i * Pi)/Pi
        W2K = 1 - (W0K + W1K);
        M2K = MT - (M0K + M1K);

        if (W2K <= 0) break;

        M2 = M2K / W2K;

        currVarB = W0K * (M0 - MT) * (M0 - MT) + W1K * (M1 - MT) * (M1 - MT) + W2K * (M2 - MT) * (M2 - MT);

        if (maxBetweenVar < currVarB)
        {
            maxBetweenVar = currVarB;
            optimalThresh1 = t1;
            optimalThresh2 = t2;
        }

    }
}

}

我其实不知道这两行是什么：

W2K = 1 - (W0K + W1K);
M2K = MT - (M0K + M1K);

以及如何将它们扩展为第三个阈值。

提前致谢。

【问题讨论】：

我看到您找到了解决方案，但这里有一些未来的信息。 Otsus 方法最适用于双峰分布。有两种非常相似的算法适用于超过 2 个集群，Expectation Maximization 和 K-means clustering 这两种方法都在 OpenCV 中。它们是 Otsus 方法的自然扩展，努力减少类方差。这两种方法甚至可以用于对 2 个或更多集群进行分组。在双峰情况下，应该产生与 otsu 方法相似的集群
@andrew，非常感谢您的评论。我一定会研究这些，看看它们是否更适合我的应用程序。我会在有更多数据后立即更新答案。
它们并不完全相同，但我认为它们会按照您想要的方式对您的直方图进行分组。就像 Otsu 的方法一样，他们假设数据是高斯的混合。 Otsu 的方法在两个类别之间找到一个阈值。而 EM 或 Kmeans 将找到实际分布（均值、方差）。使用这些属性 otsus 方法计算阈值（如您在currVarB 行中看到的那样），因此如果这是您真正想要的，您可以计算阈值，否则使用均值和方差，具体取决于您的应用程序
在http://stackoverflow.com/questions/22706742/multi-otsumulti-thresholding-with-opencv看我的评论我发现“Silva”算法有一个错误
你说得对@reexre，非常感谢。我编辑了答案。

标签： opencv image-processing computer-vision

【解决方案1】：

[编辑] 修复了@reexre 报告的错误是的，没关系。

我实现了第三个内部 for 循环来添加另一个阈值级别，并使用以下几行作为 W3K 和 M3K，尽管我仍然不知道它们为什么会这样。

W3K = 1 - (W1K + W2K);
M3K = MT - (M1K + M2K);

有效。

整个代码是，那么：

#include "Shadow01-1.cuh"

void multiThresh(double &optimalThresh1, double &optimalThresh2, double &optimalThresh3, cv::Mat &imgHist, cv::Mat &src)
{
double W0K, W1K, W2K, W3K, M0, M1, M2, M3, currVarB, maxBetweenVar, M0K, M1K, M2K, M3K, MT;
unsigned char *histogram = (unsigned char*)(imgHist.data);

int N = src.rows*src.cols;
W0K = 0;
W1K = 0;
M0K = 0;
M1K = 0;
MT = 0;
maxBetweenVar = 0;

for (int k = 0; k <= 255; k++) {
    MT += k * (histogram[k] / (double) N);
}

for (int t1 = 0; t1 <= 255; t1++)
{
    W0K += histogram[t1] / (double) N; //Pi
    M0K += t1 * (histogram[t1] / (double) N); //i * Pi
    M0 = M0K / W0K; //(i * Pi)/Pi

    W1K = 0;
    M1K = 0;

    for (int t2 = t1 + 1; t2 <= 255; t2++)
    {
        W1K += histogram[t2] / (double) N; //Pi
        M1K += t2 * (histogram[t2] / (double) N); //i * Pi
        M1 = M1K / W1K; //(i * Pi)/Pi
        W2K = 1 - (W0K + W1K);
        M2K = MT - (M0K + M1K);

        if (W2K <= 0) break;

        M2 = M2K / W2K;

        W2K = 0;
        M2K = 0;

        for (int t3 = t2 + 1; t3 <= 255; t3++)
        {
            W2K += histogram[t3] / (double) N; //Pi
            M2K += t3 * (histogram[t3] / (double) N); // i*Pi
            M2 = M2K / W2K; //(i*Pi)/Pi
            W3K = 1 - (W0K + W1K + W2K);
            M3K = MT - (M0K + M1K + M2K);

            M3 = M3K / W3K;
            currVarB = W0K * (M0 - MT) * (M0 - MT) + W1K * (M1 - MT) * (M1 - MT) + W2K * (M2 - MT) * (M2 - MT) + W3K * (M3 - MT) * (M3 - MT);

            if (maxBetweenVar < currVarB)
            {
                maxBetweenVar = currVarB;
                optimalThresh1 = t1;
                optimalThresh2 = t2;
                optimalThresh3 = t3;
            }
        }
    }
}
}

【讨论】：