用OpenCV實現Otsu算法

一、Otsu算法原理

Otsu算法（大津法或最大類間方差法）使用的是聚類的思想，把圖像的灰度數按灰度級分成2個部分，使得兩個部分之間的灰度值差異最大，每個部分之間的灰度差異最小，通過方差的計算來尋找一個合適的灰度級別來划分。 所以可以在二值化的時候采用otsu算法來自動選取閾值進行二值化。otsu算法被認為是圖像分割中閾值選取的最佳算法，計算簡單，不受圖像亮度和對比度的影響。因此，使類間方差最大的分割意味着錯分概率最小。

設t為設定的閾值。

w0	分開后前景像素點數占圖像的比例
u0	分開后前景像素點的平均灰度
w1	分開后背景像素點數占圖像的比例
u1	分開后背景像素點的平均灰度

 

 

 

 

圖像總平均灰度為： u = w0∗u0 + w1∗u1 

 

從L個灰度級遍歷 t，使得 t 為某個值的時候，前景和背景的方差最大，則 這個 t 值便是我們要求得的閾值。其中，方差的計算公式如下：

 g = wo∗(u0−u)∗(u0−u) + w1∗(u1−u)∗(u1−u) 

此公式計算量較大，可以采用：

 g = w0∗w1∗(u0−u1)∗(u0−u1) 

由於Otsu算法是對圖像的灰度級進行聚類，因此在執行Otsu算法之前，需要計算該圖像的灰度直方圖。

二、代碼實現算法

第一份代碼：來自Augusdi

#include<iostream>
#include<opencv\cv.h>
//#include<opencv2\highgui.hpp>
#include<opencv2\highgui\highgui.hpp>

int Otsu(IplImage* src);

int main()
{
    IplImage* img = cvLoadImage("E:\\house.jpg", 0);
    IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1);
    int threshold = Otsu(img);

    cvThreshold(img, dst, threshold, 255, CV_THRESH_BINARY);


    cvNamedWindow("img", 1);
    cvShowImage("img", dst);


    cvWaitKey(-1);

    cvReleaseImage(&img);
    cvReleaseImage(&dst);

    cvDestroyWindow("dst");
    return 0;
}

int Otsu(IplImage* src)
{
    int height = src->height;
    int width = src->width;
    long size = height * width;

    //histogram    
    float histogram[256] = { 0 };
    for (int m = 0; m < height; m++)
    {
        unsigned char* p = (unsigned char*)src->imageData + src->widthStep * m;
        for (int n = 0; n < width; n++)
        {
            histogram[int(*p++)]++;
        }
    }

    int threshold;
    long sum0 = 0, sum1 = 0; //存儲前景的灰度總和和背景灰度總和  
    long cnt0 = 0, cnt1 = 0; //前景的總個數和背景的總個數  
    double w0 = 0, w1 = 0; //前景和背景所占整幅圖像的比例  
    double u0 = 0, u1 = 0;  //前景和背景的平均灰度  
    double variance = 0; //最大類間方差  
    int i, j;
    double u = 0;
    double maxVariance = 0;
    for (i = 1; i < 256; i++) //一次遍歷每個像素  
    {
        sum0 = 0;
        sum1 = 0;
        cnt0 = 0;
        cnt1 = 0;
        w0 = 0;
        w1 = 0;
        for (j = 0; j < i; j++)
        {
            cnt0 += histogram[j];
            sum0 += j * histogram[j];
        }

        u0 = (double)sum0 / cnt0;
        w0 = (double)cnt0 / size;

        for (j = i; j <= 255; j++)
        {
            cnt1 += histogram[j];
            sum1 += j * histogram[j];
        }

        u1 = (double)sum1 / cnt1;
        w1 = 1 - w0; // (double)cnt1 / size;  

        u = u0 * w0 + u1 * w1; //圖像的平均灰度  
        printf("u = %f\n", u);
        //variance =  w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);  
        variance = w0 * w1 *  (u0 - u1) * (u0 - u1);
        if (variance > maxVariance)
        {
            maxVariance = variance;
            threshold = i;
        }
    }

    printf("threshold = %d\n", threshold);
    return threshold;
}

 

第二份代碼：

#include <opencv\cv.h>
#include <opencv2\highgui\highgui.hpp>

int otsu(IplImage *image)
{
    assert(NULL != image);

    int width = image->width;
    int height = image->height;
    int x = 0, y = 0;
    int pixelCount[256];
    float pixelPro[256];
    int i, j, pixelSum = width * height, threshold = 0;

    uchar* data = (uchar*)image->imageData;

    //初始化
    for (i = 0; i < 256; i++)
    {
        pixelCount[i] = 0;
        pixelPro[i] = 0;
    }

    //統計灰度級中每個像素在整幅圖像中的個數
    for (i = y; i < height; i++)
    {
        for (j = x; j <width; j++)
        {
            pixelCount[data[i * image->widthStep + j]]++;
        }
    }

    //計算每個像素在整幅圖像中的比例
    for (i = 0; i < 256; i++)
    {
        pixelPro[i] = (float)(pixelCount[i]) / (float)(pixelSum);
    }

    //經典ostu算法,得到前景和背景的分割
    //遍歷灰度級[0,255],計算出方差最大的灰度值,為最佳閾值
    float w0, w1, u0tmp, u1tmp, u0, u1, u, deltaTmp, deltaMax = 0;
    for (i = 0; i < 256; i++)
    {
        w0 = w1 = u0tmp = u1tmp = u0 = u1 = u = deltaTmp = 0;

        for (j = 0; j < 256; j++)
        {
            if (j <= i) //背景部分
            {
                //以i為閾值分類，第一類總的概率
                w0 += pixelPro[j];
                u0tmp += j * pixelPro[j];
            }
            else       //前景部分
            {
                //以i為閾值分類，第二類總的概率
                w1 += pixelPro[j];
                u1tmp += j * pixelPro[j];
            }
        }

        u0 = u0tmp / w0;        //第一類的平均灰度
        u1 = u1tmp / w1;        //第二類的平均灰度
        u = u0tmp + u1tmp;        //整幅圖像的平均灰度
        //計算類間方差
        deltaTmp = w0 * (u0 - u)*(u0 - u) + w1 * (u1 - u)*(u1 - u);
        //找出最大類間方差以及對應的閾值
        if (deltaTmp > deltaMax)
        {
            deltaMax = deltaTmp;
            threshold = i;
        }
    }
    //返回最佳閾值;
    return threshold;
}

int main(int argc, char* argv[])
{
    IplImage* srcImage = cvLoadImage("house.jpg", 0);
    assert(NULL != srcImage);

    cvNamedWindow("src");
    cvShowImage("src", srcImage);

    IplImage* biImage = cvCreateImage(cvGetSize(srcImage), 8, 1);

    //計算最佳閾值
    int threshold = otsu(srcImage);
    printf("threshold = %d\n", threshold);
    //對圖像二值化
    cvThreshold(srcImage, biImage, threshold, 255, CV_THRESH_BINARY);

    cvNamedWindow("binary");
    cvShowImage("binary", biImage);

    cvWaitKey(0);

    cvReleaseImage(&srcImage);
    cvReleaseImage(&biImage);
    cvDestroyWindow("src");
    cvDestroyWindow("binary");

    return 0;
}

 

三、代碼運行結果

代碼1運行結果：

代碼2運行結果：