首页 文章

使用openCV进行多otsu（多阈值处理）

提问于
浏览 次
9

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

关于如何进行多级阈值处理或者更确切地说是递归阈值处理的文档相当有限 . 获得原始otsu的 Value 后我该怎么办？会欣赏一些提示，我一直在玩代码，添加一个外部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;
opencv image-processing computer-vision

6 回答

  • 0

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

    multi between class variance

    请查看Huang-Yuan Huang，Lin-Wei Lin，Wu-Chih Hu，基于两阶段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） .

    original scan

    thresholded scan

    image histogram

    Otsu的原始论文：“Nobuyuki Otsu，灰度直方图的阈值选择方法，IEEE系统，人与控制论交易，1979,9：62-66”也简要提到了多阈值的扩展 . https://engineering.purdue.edu/kak/computervision/ECE661.08/OTSU_paper.pdf

    希望这可以帮助 .

    回复于
  • 2

    我已经写过一个关于otsu阈值处理如何在python中工作的例子 . 你可以在这里看到源代码：https://github.com/subokita/Sandbox/blob/master/otsu.py

    在示例中有2个变体，otsu2（）是优化版本，如维基百科页面上所见，而otsu（）是基于算法描述本身的更天真的实现 .

    如果你可以阅读python代码（在这种情况下，它们非常简单，几乎是伪代码），你可能想看看示例中的otsu（）并修改它 . 将它移植到C代码也不难 .

    回复于
  • 0

    @ Antoni4在我看来给出了最好的答案，并且非常直接地增加了关卡数量 .

    这是针对三级阈值处理：

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

        W3K = 0;
        M3K = 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 - (W1K + W2K);
            M3K = MT - (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;
            }
        }
    }
}
}
    回复于
  • 11

    @Guilherme Silva

    你的代码有一个BUG

    你必须替换：

    W3K = 0;
    M3K = 0;

    同

    W2K = 0;
    M2K = 0;

    和

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

    同

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

    ;-） 问候

    编辑（1）：[Toby Speight]我发现这个错误是通过将效果应用于不同重塑（Sizes）的相同图片，并看到输出结果彼此差异很大（甚至更改分辨率）

    W3K和M3K必须是总数减去前一个WK和MK . （我认为这与代码相似性与一个级别较少的那个）目前由于我缺乏英语我无法解释更好的方式和为什么

    说实话，我仍然不能100％确定这种方式是正确的，即使从我的输出中我也可以看出它可以提供更好的结果 . （即使有1级以上（5级灰色））你可以试试自己;-)抱歉

    我的输出：

    3 Thresholds 4 Thresholds

    回复于
  • 1

    我在这个帖子中找到了一段有用的代码 . 我正在寻找双层/浮点图像的多级Otsu实现 . 所以，我试图用双/浮点矩阵作为输入推广N级的例子 . 在我的下面的代码中，我使用armadillo库作为依赖 . 但是这个代码可以很容易地适应标准的C数组，只需用单维双重和整数数组替换vec，uvec对象，使用二维替换mat和umat . 这里使用的犰狳的另外两个函数是：vectorise和hist .

    // Input parameters:
// map - input image (double matrix)
// mask - region of interest to be thresholded
// nBins - number of bins
// nLevels - number of Otsu thresholds

#include <armadillo>
#include <algorithm>
#include <vector>

mat OtsuFilterMulti(mat map, int nBins, int nLevels) {

    mat mapr;   // output thresholded image
    mapr = zeros<mat>(map.n_rows, map.n_cols);

    unsigned int numElem = 0;
    vec threshold = zeros<vec>(nLevels);
    vec q = zeros<vec>(nLevels + 1);
    vec mu = zeros<vec>(nLevels + 1);
    vec muk = zeros<vec>(nLevels + 1);
    uvec binv = zeros<uvec>(nLevels);

    if (nLevels <= 1) return mapr;

    numElem = map.n_rows*map.n_cols;


    uvec histogram = hist(vectorise(map), nBins);

    double maxval = map.max();
    double minval = map.min();
    double odelta = (maxval - abs(minval)) / nBins;     // distance between histogram bins

    vec oval = zeros<vec>(nBins);
    double mt = 0, variance = 0.0, bestVariance = 0.0;

    for (int ii = 0; ii < nBins; ii++) {
        oval(ii) = (double)odelta*ii + (double)odelta*0.5;  // centers of histogram bins
        mt += (double)ii*((double)histogram(ii)) / (double)numElem;
    }

    for (int ii = 0; ii < nLevels; ii++) {
        binv(ii) = ii;
    }

    double sq, smuk;
    int nComb;

    nComb = nCombinations(nBins,nLevels);
    std::vector<bool> v(nBins);
    std::fill(v.begin(), v.begin() + nLevels, true);

    umat ibin = zeros<umat>(nComb, nLevels); // indices from combinations will be stored here

    int cc = 0;
    int ci = 0;
    do {
        for (int i = 0; i < nBins; ++i) {
            if(ci==nLevels) ci=0;
                if (v[i]) {
                ibin(cc,ci) = i;
                ci++;
            }
        }
        cc++;
    } while (std::prev_permutation(v.begin(), v.end()));

    uvec lastIndex = zeros<uvec>(nLevels);

    // Perform operations on pre-calculated indices
    for (int ii = 0; ii < nComb; ii++) {
        for (int jj = 0; jj < nLevels; jj++) {
            smuk = 0;
            sq = 0;
            if (lastIndex(jj) != ibin(ii, jj) || ii == 0) {
                q(jj) += double(histogram(ibin(ii, jj))) / (double)numElem;
                muk(jj) += ibin(ii, jj)*(double(histogram(ibin(ii, jj)))) / (double)numElem;
                mu(jj) = muk(jj) / q(jj);
                q(jj + 1) = 0.0;
                muk(jj + 1) = 0.0;

                if (jj>0) {
                    for (int kk = 0; kk <= jj; kk++) {
                        sq += q(kk);
                        smuk += muk(kk);
                    }
                    q(jj + 1) = 1 - sq;
                    muk(jj + 1) = mt - smuk;
                    mu(jj + 1) = muk(jj + 1) / q(jj + 1);
                }
                if (jj>0 && jj<(nLevels - 1)) {
                    q(jj + 1) = 0.0;
                    muk(jj + 1) = 0.0;
                }

                lastIndex(jj) = ibin(ii, jj);
            }
        }

        variance = 0.0;
        for (int jj = 0; jj <= nLevels; jj++) {
            variance += q(jj)*(mu(jj) - mt)*(mu(jj) - mt);
        }    

        if (variance > bestVariance) {
            bestVariance = variance;
            for (int jj = 0; jj<nLevels; jj++) {
                threshold(jj) = oval(ibin(ii, jj));
            }
        }    
    }

    cout << "Optimized thresholds: ";
    for (int jj = 0; jj<nLevels; jj++) {
        cout << threshold(jj) << " ";
    }
    cout << endl;

    for (unsigned int jj = 0; jj<map.n_rows; jj++) {
        for (unsigned int kk = 0; kk<map.n_cols; kk++) {
            for (int ll = 0; ll<nLevels; ll++) {
                if (map(jj, kk) >= threshold(ll)) {
                    mapr(jj, kk) = ll+1;
                }
            }
        }
    }

    return mapr;
}


int nCombinations(int n, int r) {

    if (r>n) return 0;
    if (r*2 > n) r = n-r;
    if (r == 0) return 1;

    int ret = n;
    for( int i = 2; i <= r; ++i ) {
        ret *= (n-i+1);
        ret /= i;
    }
    return ret;
}
    回复于
  • 1

    以下是 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)

    有关完整论文和更多信息，请访问this link .

    回复于

相关问题