图像分割：Otsu大津算法阈值选择
1、Otsu算法原理：
Otsu算法使用的是聚类的思想：
①把图像的灰度数按灰度级分成2个部分，使得两个部分之间的灰度值差异最大，每个部分之间的灰度差异最小；
②通过方差的计算来寻找一个合适的灰度级别来划分；
特点：
它是按图像的灰度特性,将图像分成背景和前景两部分；
因方差是灰度分布均匀性的一种度量,背景和前景之间的类间方差越大,说明构成图像的两部分的差别越大,当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小；
因此,使类间方差最大的分割意味着错分概率最小。
2、对图像I(x,y)，前景和背景的分割阈值记作T，
前景像素点数占整幅图像的比例为ω0，其平均灰度μ0；
背景像素点数占整幅图像的比例为ω1，其平均灰度μ1；
图像的总平均灰度记为μ=ω0∗μ0 +ω1∗μ1；
类间方差记为g。
假设：
背景较暗，且图像的大小为M×N，
图像中，
像素灰度值小于阈值T的像素个数记作N0，
像素灰度值大于阈值T的像素个数记作N1，
则有：
ω0 = N0 / M×N ； (1)
ω1 = N1 / M×N ； (2)
N0 + N1 = M×N ； (3)
ω0+ω1=1； (4)
μ=ω0*μ0+ω1*μ1； (5)
g =ω0(μ0-μ)^2+ω1(μ1-μ)^2；(6)
将式(5)代入式(6),得到等价公式:
g=ω0*ω1*(μ0-μ1)^2; (7) 这就是类间方差
采用遍历的方法得到使类间方差g最大的阈值T,即为所求。
【注】：由于Otsu算法是对图像的灰度级进行聚类，因此在执行Otsu算法之前，需要计算该图像的灰度直方图；
3、matlab算法实现：
matlab中函数graythresh()即使用大津法求得分割阈值Ｔ；
用法如下：
T = graythresh(img);%求阈值
BW = im2bw(img,T);%二值化
大津法的形象理解：
对于直方图有两个峰值的图像，大津法求得的Ｔ近似等于两个峰值之间的低谷；
imhist(img);
.........................T = graythresh(img);
下图为图像的直方图，使大津法求得的Ｔ＝0.5294，转换在［0,255］之间为134.9970，正好是两个峰值之间低谷的位置；
4、OpenCV算法实现：
根据大津法思路编写程序：
代码一：
#include <opencv2\opencv.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\features2d\features2d.hpp>
#include <opencv2\core\core.hpp>
using namespace std;
using namespace cv;
int OtsuThreshold(IplImage *frame);
void main()
{
IplImage *raw_img = cvLoadImage("4-1.jpg",0);
cvShowImage("0",raw_img);
double t = (double)getTickCount();
int m_threshold = OtsuThreshold(raw_img);
//返回该处代码执行所耗的时间，单位为秒
t = ((double)getTickCount() - t)/getTickFrequency();
cout<<t*1000<<"ms"<<endl;
for (int j=0; j<raw_img->width;j++)//lie
{
for (int i=0;i<raw_img->height;i++)
{
unsigned char *p_data = (unsigned char*)raw_img->imageData + i*raw_img->widthStep;
if (p_data[j] < m_threshold)
{
p_data[j] = 0;
}
else
p_data[j]=255;
}
}
cvShowImage("1",raw_img);
cvWaitKey(0);
}
int OtsuThreshold(IplImage *frame)
{
int width = frame->width;//lie
int height = frame->height;//hang
const int GrayScale = 256;
int pixelCount[GrayScale];
float pixelPro[GrayScale];
int i;
for (i = 0; i < GrayScale; i++)//数组初始化
{
pixelCount[i] = 0;
pixelPro[i] = 0;
}
int j, pixelSum = width * height, threshold = 0;
unsigned char* data = (unsigned char*)frame->imageData; //指向像素数据的指针
for (i = 0; i < height; i++)//统计灰度级中每个像素在整幅图像中的个数
{
for (j = 0; j < width; j++)
{
pixelCount[(int)data[i * width + j]]++; //将像素值作为计数数组的下标
}
}
//计算每个像素在整幅图像中的比例
float maxPro = 0.0;
int kk = 0;
for (i = 0; i < GrayScale; i++)
{
pixelPro[i] = (float)pixelCount[i] / pixelSum;
if (pixelPro[i] > maxPro)
{
maxPro = pixelPro[i];//最大比例
kk = i;
}
}
//遍历灰度级[0,255]
float w0, w1, u0tmp, u1tmp, u0, u1, u, deltaTmp, deltaMax = 0;
for (i = 0; i < GrayScale; i++) //遍历灰度等级，i为阈值
{
w0 = w1 = u0tmp = u1tmp = u0 = u1 = u = deltaTmp = 0;
for (j = 0; j < GrayScale; j++)
{
if (j <= i) //背景部分
{
w0 += pixelPro[j];//背景像素点占整个图像的比例
u0tmp += j * pixelPro[j];
}
else //前景部分
{
w1 += pixelPro[j];//前景像素点占整个图像的比例
u1tmp += j * pixelPro[j];
}
}
u0 = u0tmp / w0;//背景平均灰度μ0
u1 = u1tmp / w1;//前景平均灰度μ1
u = u0tmp + u1tmp;
deltaTmp = w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);//g =ω0(μ0-μ)^2+ω1(μ1-μ)^2;类间方差
if (deltaTmp > deltaMax)
{
deltaMax = deltaTmp;
threshold = i;
}
}
return threshold;
}
5、处理一张图片花费时间多，非常低效，但是是一种非常直接的，显而易见的实现方法；
代码二：
int otsu(IplImage* src)
{
int Width = src->width;
int Height = src->height;
CvScalar cs = cvSum(src); //图像灰度值累加和
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++) //lie
{
for(int j=0; j < Height;j++) //hang
{
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;
}
6、代码三：
借用opencv中的直方图函数，大大减少循环数；
【注】：
预先宏定义：
#define cvQueryHistValue_1D( hist, idx0 ) ((float)cvGetReal1D( (hist)->bins, (idx0)))
int otsu_1(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;
}
7、代码四：
一次循环实现，其中的公式需要自己用笔推导一下，非常精妙！
【注】：
预先宏定义：
#define cvQueryHistValue_1D( hist, idx0 ) ((float)cvGetReal1D( (hist)->bins, (idx0)))
int otsu_2(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;
}
