用OpenCV实现Otsu算法

正文

最近在学习图像分割反面的知识，在冈萨雷斯的那本书上看到Otsu算法，身边的同学都是用Matlab来实现这个算法。我觉得Matlab写得话，但是代码的效率应该不会高。于是又恶补了一些OpenCV的一些基本知识，然后看了Augusdi的博客，分析了一下他的代码并附上，第二份代码是来自某一位大牛的代码，写得更加易懂。

---

算法的介绍

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算法是对图像的灰度级进行聚类，so 在执行otsu算法之前，需要计算该图像的灰度直方图。

---

代码实现算法

第一份代码：来自Augusdi

#include<iostream>#include<opencv\cv.h>#include<opencv2\highgui.hpp>int Otsu(IplImage* src);int main(){    IplImage* img = cvLoadImage("E:\\pp.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("E:\\pp.jpg", 0);    assert(NULL != srcImage);    cvNamedWindow("src");    cvShowImage("src", srcImage);    IplImage* biImage = cvCreateImage(cvGetSize(srcImage), 8, 1);    //计算最佳阈值    int threshold = otsu(srcImage);    //对图像二值化    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;}

---

代码运行结果