使用openCV进行多otsu(多阈值处理)

Question

我正在尝试用otsu进行多阈值处理.我目前使用的方法实际上是通过最大化类间差异,我已经设法获得与OpenCV库相同的阈值.但是,这只是通过运行otsu方法一次.

关于如何进行多级阈值处理或者更确切地说是递归阈值处理的文档相当有限.获得原始otsu的价值后我该怎么办？会欣赏一些提示,我一直在玩代码,添加一个外部for循环,但计算的下一个值对于任何给定的图像总是254 :(

我的代码如果需要:

//compute histogram first
cv::Mat imageh; //image edited to grayscale for histogram purpose
//imageh=image; //to delete and uncomment below;
cv::cvtColor(image, imageh, CV_BGR2GRAY);

int histSize[1] = {256}; // number of bins
float hranges[2] = {0.0, 256.0}; // min andax pixel value
const float* ranges[1] = {hranges};
int channels[1] = {0}; // only 1 channel used

cv::MatND hist;
// Compute histogram
calcHist(&imageh, 1, channels, cv::Mat(), hist, 1, histSize, ranges);

IplImage* im = new IplImage(imageh);//assign the image to an IplImage pointer
IplImage* finalIm = cvCreateImage(cvSize(im->width, im->height), IPL_DEPTH_8U, 1);
double otsuThreshold= cvThreshold(im, finalIm, 0, 255, cv::THRESH_BINARY | cv::THRESH_OTSU );

cout<<"opencv otsu gives "<<otsuThreshold<<endl;

int totalNumberOfPixels= imageh.total();
cout<<"total number of Pixels is " <<totalNumberOfPixels<< endl;


float sum = 0;
for (int t=0 ; t<256 ; t++) 
{
    sum += t * hist.at<float>(t);
}
cout<<"sum is "<<sum<<endl;

float sumB = 0; //sum of background
int wB = 0; // weight of background
int wF = 0; //weight of foreground

float varMax = 0;
int threshold = 0;

//run an iteration to find the maximum value of the between class variance(as between class variance shld be maximise)
for (int t=0 ; t<256 ; t++) 
{
       wB += hist.at<float>(t);               // Weight Background
       if (wB == 0) continue;

       wF = totalNumberOfPixels - wB;                 // Weight Foreground
       if (wF == 0) break;

       sumB += (float) (t * hist.at<float>(t));

       float mB = sumB / wB;            // Mean Background
       float mF = (sum - sumB) / wF;    // Mean Foreground

       // Calculate Between Class Variance
       float varBetween = (float)wB * (float)wF * (mB - mF) * (mB - mF);

       // Check if new maximum found
       if (varBetween > varMax) {
          varMax = varBetween;
          threshold = t;
       }
}

       cout<<"threshold value is: "<<threshold;

Answer 1

为了将Otsu的阈值方法扩展到多级阈值处理,类间方差方程变为:

请查看黄登元,林大伟,胡武胡,基于两阶段Otsu方法的自动多级阈值,通过Valley Estimation进行聚类确定,Int.Journal of Innovative Computing,2011,7:5631-5644获取更多信息.

http://www.ijicic.org/ijicic-10-05033.pdf

这是我的Otsu Multi的C#实现2阈值:

/* Otsu (1979) - multi */

Tuple < int, int > otsuMulti(object sender, EventArgs e) {
    //image histogram
    int[] histogram = new int[256];

    //total number of pixels
    int N = 0;

    //accumulate image histogram and total number of pixels
    foreach(int intensity in image.Data) {
        if (intensity != 0) {
            histogram[intensity] += 1;
            N++;
        }
    }

    double W0K, W1K, W2K, M0, M1, M2, currVarB, optimalThresh1, optimalThresh2, maxBetweenVar, M0K, M1K, M2K, MT;

    optimalThresh1 = 0;
    optimalThresh2 = 0;

    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;
            }
        }
    }

    return new Tuple(optimalThresh1, optimalThresh2);
}

这是我用以上代码对土壤图像扫描进行阈值处理得到的结果:

(T1 = 110,T2 = 147).

Otsu的原始论文:"Nobuyuki Otsu,灰度直方图的阈值选择方法,IEEE系统,人类和控制论交易,1979,9:62-66"也简要提到了多重阈值的扩展.

https://engineering.purdue.edu/kak/computervision/ECE661.08/OTSU_paper.pdf

希望这可以帮助.

Answer 2

这是python (>3.0)中“n”个阈值的简单通用方法：

# developed by- SUJOY KUMAR GOSWAMI
# source paper- https://people.ece.cornell.edu/acharya/papers/mlt_thr_img.pdf

import cv2
import numpy as np
import math

img = cv2.imread('path-to-image')
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
a = 0
b = 255
n = 6 # number of thresholds (better choose even value)
k = 0.7 # free variable to take any positive value
T = [] # list which will contain 'n' thresholds

def sujoy(img, a, b):
    if a>b:
        s=-1
        m=-1
        return m,s

    img = np.array(img)
    t1 = (img>=a)
    t2 = (img<=b)
    X = np.multiply(t1,t2)
    Y = np.multiply(img,X)
    s = np.sum(X)
    m = np.sum(Y)/s
    return m,s

for i in range(int(n/2-1)):
    img = np.array(img)
    t1 = (img>=a)
    t2 = (img<=b)
    X = np.multiply(t1,t2)
    Y = np.multiply(img,X)
    mu = np.sum(Y)/np.sum(X)

    Z = Y - mu
    Z = np.multiply(Z,X)
    W = np.multiply(Z,Z)
    sigma = math.sqrt(np.sum(W)/np.sum(X))

    T1 = mu - k*sigma
    T2 = mu + k*sigma

    x, y = sujoy(img, a, T1)
    w, z = sujoy(img, T2, b)

    T.append(x)
    T.append(w)

    a = T1+1
    b = T2-1
    k = k*(i+1)

T1 = mu
T2 = mu+1
x, y = sujoy(img, a, T1)
w, z = sujoy(img, T2, b)    
T.append(x)
T.append(w)
T.sort()
print(T)

如需全文和更多信息，请访问此链接。

Answer 3

我之前写过一个关于 otsu 阈值如何在 python 中工作的示例。您可以在这里查看源代码：https ://github.com/subokita/Sandbox/blob/master/otsu.py

在示例中，有 2 个变体，otsu2() 是优化版本，如维基百科页面所示，otsu() 是基于算法描述本身的更简单的实现。

如果您可以阅读 python 代码（在本例中，它们非常简单，几乎是伪代码），您可能需要查看示例中的 otsu() 并修改它。将其移植到 C++ 代码也不难。