OTSU算法:就是计算出灰度图最佳阈值的算法
1.先对灰度图进行直方图计算并归一化处理,得到0-255之间每个像素在灰度图中出现的概率,即表示为某个像素在灰度图中出现了n个,灰度图总的像素点为N个,则这个像素的出现概率为Pi=n/N
2.每个灰度图可以由阈值k将灰度图分为A,B两大类,很容易得到A,B类在灰度图中的出现概率以及灰度均值
3.计算灰度图A,B类得类间方差,在最佳阈值K处,求得的类间方差最大,也就是类间方差最大的那个时刻的阈值就为灰度图的最佳阈值
代码实现:
以下三种算法时间复杂度依次降低,前两种算法自己所写,最后一种算法从网上看到大神写的,膜拜!当然如果你想偷懒,可是更直接点用cvThreshold一步实现cvThreshold(srcCopy,thresholdImage,0,255,CV_THRESH_OTSU),花费时间与算法三相当。
算法一:
处理一张640*480的图片花费时间大概在8000毫秒左右,非常低效,但是是一种非常直接的,显而易见的实现方法。
int otsu(IplImage* src)
{
CvScalar cs = cvSum(src); //图像总灰度值
int Width = src->width;
int Height = src->height;
double U_t = 1.0*cs.val[0]/(Width*Height); //图像的平均灰度
int threshold = 0;
double delta = 0;
for(int k=0; k<256; k++)
{
double sumOfGray = 0;
int num = 0;
for(int i=0; i < Width; i++)
{
for(int j=0; j < Height;j++)
{
int t = cvGet2D(src, j,i).val[0];
if(t < k)
{
num++; //统计灰度小于k的像素点的个数
sumOfGray += t; //灰度值小于k的像素点的灰度和
}
}
}
double w0 = 1.0*num/(Width*Height); //灰度值小于阈值k的像素的概率
double w1 = 1.0 - w0; //灰度值大于k像素的概率
double u0 = 1.0*sumOfGray/num; //灰度值小于阈值k的像素的平均灰度值
double u1 = 1.0*(cs.val[0]-sumOfGray)/(Width*Height-num); //灰度值大于阈值k的像素的平均灰度值
double delta_tmp = w0*pow(u0-U_t,2) + w1*pow(u1-U_t,2); //类间方差
if(delta_tmp > delta)
{
delta = delta_tmp;
threshold = k;
}
}
cout<<cs.val[0]<<endl;
cout<<U_t<<endl;
return threshold;
}
算法二:
借用opencv中的直方图函数,大大减少循环数,同一张图片花费时间大概在16毫秒左右。
int otsu(IplImage* src)
{
int hist_size = 256; //直方图尺寸
int hist_height = 256;
float range[] = {0,255}; //灰度级的范围
float* ranges[]={range};
//创建一维直方图,统计图像在[0 255]像素的均匀分布
CvHistogram* gray_hist = cvCreateHist(1,&hist_size,CV_HIST_ARRAY,ranges,1);
//计算灰度图像的一维直方图
cvCalcHist(&src,gray_hist,0,0);
//归一化直方图
cvNormalizeHist(gray_hist,1.0);
int Width = src->width;
int Height = src->height;
int threshold = 0;
double delta = 0;
CvScalar sumOfGray = cvSum(src);
double U_t = sumOfGray.val[0]/(Width*Height);
for(int k=0; k<256; k++)
{
double u0 = 0, u1 = 0, w0 = 0, w1 = 0, delta_tmp = 0, u00 = 0, u11 = 0;
for(int i=0; i<k; i++)
{
u00 += cvQueryHistValue_1D(gray_hist,i)*i; //灰度小于阈值k的像素的平均灰度值
w0 += cvQueryHistValue_1D(gray_hist,i); //灰度小于阈值k的像素的概率
}
u0 = u00/w0;
for(int j=k; j<256; j++)
{
u11 += cvQueryHistValue_1D(gray_hist,j)*j; //灰度大于阈值k的像素的平均灰度值
w1 += cvQueryHistValue_1D(gray_hist,j); //灰度小于阈值k的像素的概率
}
u1 = u11/w1;
delta_tmp = w0*pow(u0-U_t,2) + w1*pow(u1-U_t,2); //类间方差
if(delta_tmp > delta)
{
delta = delta_tmp;
threshold = k;
}
}
cvReleaseHist(&gray_hist);
return threshold;
}
算法三:
一次循环实现,其中的公式需要自己用笔推导一下,非常精妙!完成一次运算的时间毫秒已经不能计量了,显示为0!
int otsu(IplImage* src)
{
int hist_size = 256; //直方图尺寸
int hist_height = 256;
float range[] = {0,255}; //灰度级的范围
float* ranges[]={range};
//创建一维直方图,统计图像在[0 255]像素的均匀分布
CvHistogram* gray_hist = cvCreateHist(1,&hist_size,CV_HIST_ARRAY,ranges,1);
//计算灰度图像的一维直方图
cvCalcHist(&src,gray_hist,0,0);
//归一化直方图
cvNormalizeHist(gray_hist,1.0);
int Width = src->width;
int Height = src->height;
int threshold = 0;
double delta = 0;
double U_t = 0;
for(int m=0; m<256; m++)
{
U_t += cvQueryHistValue_1D(gray_hist,m)*m;
}
double u = 0, w = 0;
for(int k=0; k<256; k++)
{
u += cvQueryHistValue_1D(gray_hist,k)*k; //
w += cvQueryHistValue_1D(gray_hist,k); //灰度大于阈值k的像素的概率
double t = U_t * w - u;
double delta_tmp = t * t / (w * (1 - w) );
if(delta_tmp > delta)
{
delta = delta_tmp;
threshold = k;
}
}
cvReleaseHist(&gray_hist);
return threshold;
}