opencv实现c++的otsu自适应阈值分割的算法描述

otsu算法选择使类间方差最大的灰度值为阈值，具有很好的效果
算法具体描述见otsu论文，或冈萨雷斯著名的数字图像处理那本书
这里给出程序流程：
1、计算直方图并归一化histogram
2、计算图像灰度均值avgValue.
3、计算直方图的零阶w[i]和一级矩u[i]
4、计算并找到最大的类间方差（between-class variance）
variance[i]=(avgValue*w[i]-u[i])*(avgValue*w[i]-u[i])/(w[i]*(1-w[i]))
对应此最大方差的灰度值即为要找的阈值
5、用找到的阈值二值化图像

我在代码中做了一些优化，所以算法描述的某些地方跟程序并不一致

otsu代码，先找阈值，继而二值化
// implementation of otsu algorithm

// author: onezeros(@yahoo.cn)

// reference: Rafael C. Gonzalez. DigitalImage Processing Using MATLAB

void cvThresholdOtsu(IplImage* src,IplImage* dst)

{

   int height=src->height;

   int width=src->width;

 

   //histogram

   float histogram[256]= {0};

   for(int i=0; i<height; i++)

    {

       unsigned char* p=(unsigned char*)src->imageData+src->widthStep*i;

       for(int j=0; j<width; j++)

       {

           histogram[*p++]++;

       }

    }

   //normalize histogram

   int size=height*width;

   for(int i=0; i<256; i++)

    {

       histogram[i]=histogram[i]/size;

    }

 

   //average pixel value

   float avgValue=0;

   for(int i=0; i<256; i++)

    {

       avgValue+=i*histogram[i];

    }

 

   int threshold;

   float maxVariance=0;

   float w=0,u=0;

   for(int i=0; i<256; i++)

    {

       w+=histogram[i];

       u+=i*histogram[i];

 

       float t=avgValue*w-u;

       float variance=t*t/(w*(1-w));

       if(variance>maxVariance)

       {

           maxVariance=variance;

           threshold=i;

       }

    }

 

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

}

 

更多情况下我们并不需要对每一帧都是用otsu寻找阈值，于是可以先找到阈值，然后用找到的阈值处理后面的图像。下面这个函数重载了上面的，返回值就是阈值。只做了一点改变

// implementation of otsu algorithm

// author: onezeros(@yahoo.cn)

// reference: Rafael C. Gonzalez. DigitalImage Processing Using MATLAB

int cvThresholdOtsu(IplImage* src)

{

   int height=src->height;

   int width=src->width;

 

   //histogram

   float histogram[256]= {0};

   for(int i=0; i<height; i++)

    {

       unsigned char* p=(unsigned char*)src->imageData+src->widthStep*i;

       for(int j=0; j<width; j++)

       {

           histogram[*p++]++;

       }

    }

   //normalize histogram

   int size=height*width;

   for(int i=0; i<256; i++)

    {

       histogram[i]=histogram[i]/size;

    }

 

   //average pixel value

   float avgValue=0;

   for(int i=0; i<256; i++)

    {

       avgValue+=i*histogram[i];

    }

 

   int threshold;

   float maxVariance=0;

   float w=0,u=0;

   for(int i=0; i<256; i++)

    {

       w+=histogram[i];

       u+=i*histogram[i];

 

       float t=avgValue*w-u;

       float variance=t*t/(w*(1-w));

       if(variance>maxVariance)

       {

           maxVariance=variance;

           threshold=i;

       }

    }

 

   return threshold;

}

 

 

我在手的自动检测中使用这个方法，效果很好。

下面是使用上述两个函数的简单的主程序，可以试运行一下，如果处理视频，要保证第一帧时，手要在图像中。

 

#include <cv.h>

#include <cxcore.h>

#include <highgui.h>

#pragma comment(lib,"cv210d.lib")

#pragmacomment(lib,"cxcore210d.lib")

#pragmacomment(lib,"highgui210d.lib")

 

#include <iostream>

using namespace std;

 

int main(int argc, char** argv)

{

#ifdef VIDEO //video process 

   CvCapture* capture=cvCreateCameraCapture(-1);

   if (!capture)

    {

       cout<<"failed to open camera"<<endl;

       exit(0);

    }

 

   int threshold=-1;

   IplImage* img;

   while (img=cvQueryFrame(capture))

    {

       cvShowImage("video",img);

       cvCvtColor(img,img,CV_RGB2YCrCb);

 

       IplImage* imgCb=cvCreateImage(cvGetSize(img),8,1);

       cvSplit(img,NULL,NULL,imgCb,NULL);

       if (threshold<0)

       {

           threshold=cvThresholdOtsu(imgCb);

       }

       //cvThresholdOtsu(imgCb,imgCb);

       cvThreshold(imgCb,imgCb,threshold,255,CV_THRESH_BINARY);

       cvErode(imgCb,imgCb);

       cvDilate(imgCb,imgCb);

 

       cvShowImage("object",imgCb);

       cvReleaseImage(&imgCb);

 

       if (cvWaitKey(3)==27) //esc

       {

           break;

       }

    }

 

   cvReleaseCapture(&capture);

 

#else //single image process 

   const char* filename=(argc>=2?argv[1]:"cr.jpg");

   IplImage* img=cvLoadImage(filename,CV_LOAD_IMAGE_GRAYSCALE);

 

   cvThresholdOtsu(img,img);

   cvShowImage( "src", img );

   char buf[256];

   sprintf_s(buf,256,"%s.otsu.jpg",filename);

   cvSaveImage(buf,img);

 

   cvErode(img,img);

   cvDilate(img,img);

   cvShowImage( "dst", img );

   sprintf_s(buf,256,"%s.otsu.processed.jpg",filename);

   cvSaveImage(buf,img);

 

   cvWaitKey(0);

#endif

 

   return 0;

}

#include <cv.h>

#include <cxcore.h>

#include <highgui.h>

#pragma comment(lib,"cv210d.lib")

#pragmacomment(lib,"cxcore210d.lib")

#pragmacomment(lib,"highgui210d.lib")

#include <iostream>

using namespace std;

int main(int argc, char** argv)

{

#ifdef VIDEO //video process

   CvCapture* capture=cvCreateCameraCapture(-1);

   if (!capture)

    {

       cout<<"failed to open camera"<<endl;

       exit(0);

    }

   int threshold=-1;

   IplImage* img;

   while (img=cvQueryFrame(capture))

    {

       cvShowImage("video",img);

       cvCvtColor(img,img,CV_RGB2YCrCb);

       IplImage* imgCb=cvCreateImage(cvGetSize(img),8,1);

       cvSplit(img,NULL,NULL,imgCb,NULL);

       if (threshold<0)

       {

           threshold=cvThresholdOtsu(imgCb);

       }

       //cvThresholdOtsu(imgCb,imgCb);

       cvThreshold(imgCb,imgCb,threshold,255,CV_THRESH_BINARY);

       cvErode(imgCb,imgCb);

       cvDilate(imgCb,imgCb);

 

       cvShowImage("object",imgCb);

       cvReleaseImage(&imgCb);

       if (cvWaitKey(3)==27) //esc

       {

           break;

       }

    }

   cvReleaseCapture(&capture);

 

#else //single image process

   const char* filename=(argc>=2?argv[1]:"cr.jpg");

   IplImage* img=cvLoadImage(filename,CV_LOAD_IMAGE_GRAYSCALE);

   cvThresholdOtsu(img,img);

   cvShowImage( "src", img );

   char buf[256];

   sprintf_s(buf,256,"%s.otsu.jpg",filename);

   cvSaveImage(buf,img);

   cvErode(img,img);

   cvDilate(img,img);

   cvShowImage( "dst", img );

   sprintf_s(buf,256,"%s.otsu.processed.jpg",filename);

   cvSaveImage(buf,img);

   cvWaitKey(0);

#endif

 

   return 0;

}