otsu自适应阈值分割的算法描述和opencv实现,及其在肤色检测中的应用

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(@https://www.360docs.net/doc/007074646.html,)

// reference: Rafael C. Gonzalez. Digital Image 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

{

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

for(int j=0; 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);

}

// implementation of otsu algorithm

// author: onezeros(@https://www.360docs.net/doc/007074646.html,)

// reference: Rafael C. Gonzalez. Digital Image 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

{

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

for(int j=0; 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(@https://www.360docs.net/doc/007074646.html,)

// reference: Rafael C. Gonzalez. Digital Image 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

{

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

for(int j=0; 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;

}

// implementation of otsu algorithm

// author: onezeros(@https://www.360docs.net/doc/007074646.html,)

// reference: Rafael C. Gonzalez. Digital Image 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

{

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

for(int j=0; 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

#include

#include

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

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

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

#include

using namespace std;

int main(int argc, char** argv)

{

#ifdef VIDEO //video process

CvCapture* capture=cvCreateCameraCapture(-1);

if (!capture)

{

cout<<"failed to open camera"<

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

#include

#include

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

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

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

#include

using namespace std;

int main(int argc, char** argv)

{

#ifdef VIDEO //video process

CvCapture* capture=cvCreateCameraCapture(-1);

if (!capture)

{

cout<<"failed to open camera"<

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;

}