基于阈值的分割算法
总结一下基于阈值的图像分割的几种算法,如下:
1. 双峰法阈值分割
Prewitt 等人于六十年代中期提出的直方图双峰法(也称 mode 法) 是典型的全局单阈值分割方法。该方法的基本思想是:假设图像中有明显的目标和背景,则其灰度直方图呈双峰分布,当灰度级直方图具有双峰特性时,选取两峰之间的谷对应的灰度级作为阈值。如果背景的灰度值在整个图像中可以合理地看作为恒定,而且所有物体与背景都具有几乎相同的对比度,那么,选择一个正确的、固定的全局阈值会有较好的效果。
int Thresh::GetMinimumThreshold(vector<int>&HistGram)
{
int Y, Iter = 0;
vector<int> HistGramC(256); // 基于精度问题,一定要用浮点数来处理,否则得不到正确的结果
vector<int> HistGramCC(256); // 求均值的过程会破坏前面的数据,因此需要两份数据
for (Y = 0; Y < 256; Y++)
{
HistGramC[Y] = HistGram[Y];
HistGramCC[Y] = HistGram[Y];
}
// 通过三点求均值来平滑直方图
while (IsDimodal(HistGramCC) == false) // 判断是否已经是双峰的图像了
{
HistGramCC[0] = (HistGramC[0] + HistGramC[0] + HistGramC[1]) / 3; // 第一点
for (Y = 1; Y < 255; Y++)
HistGramCC[Y] = (HistGramC[Y - 1] + HistGramC[Y] + HistGramC[Y + 1]) / 3; // 中间的点
HistGramCC[255] = (HistGramC[254] + HistGramC[255] + HistGramC[255]) / 3; // 最后一点
HistGramC = HistGramCC;
Iter++;
if (Iter >= 1000) return -1; // 直方图无法平滑为双峰的,返回错误代码
}
// 阈值极为两峰之间的最小值
bool Peakfound = false;
for (Y = 1; Y < 255; Y++)
{
if (HistGramCC[Y - 1] < HistGramCC[Y] && HistGramCC[Y + 1] < HistGramCC[Y]) Peakfound = true;
if (Peakfound == true && HistGramCC[Y - 1] >= HistGramCC[Y] && HistGramCC[Y + 1] >= HistGramCC[Y])
return Y - 1;
}
return -1;
}
其中IsDimodal函数为判断直方图是否是双峰的函数,代码如下:
bool IsDimodal(vector<int>&HistGram)
{
//对直方图的峰进行计数,只有在峰数位2才为双峰
int Count = 0;
for (int i = 1; i < 255; i++)
{
if (HistGram[i - 1] < HistGram[i] && HistGram[i + 1] < HistGram[i])
{
Count++;
if (Count > 2)
return false;
}
}
if (Count == 2)
return true;
else
return false;
}
直方图计算代码如下:
void Calcu_HistGram(Mat&img,vector<int>&HistGram)
{
if (img.channels() == 3)
cvtColor(img, img, COLOR_BGR2GRAY);
img.convertTo(img, CV_8UC1);
for (int k = 0; k < 256; k++)
{
HistGram[k] = 0;
}
int sum = 0;
for (int i = 0; i < img.rows; ++i)
{
for (int j = 0; j < img.cols; ++j)
{
HistGram[img.at<uchar>(i, j)]++;
sum++;
}
}
//cout << sum << endl;
}
2.迭代阈值分割
算法原理如下:
int GetIterativeBestThreshold(vector<int>&HistGram HistGram)
{
int X, Iter = 0;
int MeanValueOne, MeanValueTwo, SumOne, SumTwo, SumIntegralOne, SumIntegralTwo;
int MinValue, MaxValue;
int Threshold, NewThreshold;
for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;
for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;
if (MaxValue == MinValue) return MaxValue; // 图像中只有一个颜色
if (MinValue + 1 == MaxValue) return MinValue; // 图像中只有二个颜色
Threshold = MinValue;
NewThreshold = (MaxValue + MinValue) >> 1;
while (Threshold != NewThreshold) // 当前后两次迭代的获得阈值相同时,结束迭代
{
SumOne = 0; SumIntegralOne = 0;
SumTwo = 0; SumIntegralTwo = 0;
Threshold = NewThreshold;
//根据阈值将图像分割成目标和背景两部分,求出两部分的平均灰度值
for (X = MinValue; X <= Threshold; X++)
{
SumIntegralOne += HistGram[X] * X;
SumOne += HistGram[X];
}
MeanValueOne = SumIntegralOne / SumOne;
for (X = Threshold + 1; X <= MaxValue; X++)
{
SumIntegralTwo += HistGram[X] * X;
SumTwo += HistGram[X];
}
MeanValueTwo = SumIntegralTwo / SumTwo;
NewThreshold = (MeanValueOne + MeanValueTwo) >> 1; //求出新的阈值
Iter++;
if (Iter >= 1000) return -1;
}
return Threshold;
}
3.自适应阈值分割
在许多情况下,物体和背景的对比度在图象中不是各处一样的,这时很难用统一的一个阈值将物体与背景分开。这时可以根据图象的局部特征分别采用不同的阈值进行分割。实际处理时,需要按照具体问题将图象分成若干子区域分别选择阈值,或者动态地根据一定的邻域范围选择每点处的阈值,进行图象分割。
大津法(OTSU),最大类间方差法是由日本学者大津于1979年提出的,是一种自适应的阈值确定的方法,又叫大津法,简称OTSU。它是按图像的灰度特性,将图像分成背景和目标2部分。背景和目标之间的类间方差越大,说明构成图像的2部分的差别越大,当部分目标错分为背景或部分背景错分为目标都会导致2部
分差别变小。因此,使类间方差最大的分割意味着错分概率最小。
对于图像I(x,y),前景(即目标)和背景的分割阈值记作T, 属于前景的像素点数占整幅图像的比例记为ω0,其平均灰度μ0;背景像素点数占整幅图像的比例为ω1,其平均灰度为μ1。图像的总平均灰度记为μ,类间方差记为g。
假设图像的背景较暗,并且图像的大小为M×N,
图像中像素的灰度值小于阈值T的像素个数记作N0,像素灰度大于阈值T的像素个数记作N1,则有:
w0=N0/MxN (1)
w1=N1/MxN (2)
N0+N1=MxN (3)
w0+w1=1 (4)
u=w0*u0+w1*u1 (5)
g=w0(u0-u)^2+w1(u1-u)^2 (6)
将式(5)代入式(6),得到等价公式:
g=w0*w1*(u0-u1)^2 (7)
采用遍历的方法得到使类间方差最大的阈值T,即为所求。
/*
parameter: *image --- buffer for image
rows, cols --- size of image
x0, y0, dx, dy --- region of vector used for computing threshold
vvv --- debug option, is 0, no debug information outputed
*/
/*======================================================================*/
/* OTSU global thresholding routine */
/* takes a 2D unsigned char array pointer, number of rows, and */
/* number of cols in the array. returns the value of the threshold */
/*======================================================================*/
int Otsu_Segment(Mat&src, Mat&dst)
{
if(src.channels() == 3)
cvtColor(src, src, COLOR_BGR2GRAY);
Mat dst = src.clone();
const int GrayScale = 256;
int width = src.cols;
int height = src.rows;
int pixelCount[GrayScale];
float pixelPro[GrayScale];
int pixelSum = width * height, threshold = 0;
uchar*data = src.data;
for (int i = 0; i < GrayScale; i++)
{
pixelCount[i] = 0;
pixelPro[i] = 0;
}
//统计灰度级中每个像素在整幅图像中的个数
for (int i = 0; i < height; ++i)
{
for (int j = 0; j < width; ++j)
{
pixelCount[data[i*width + j]]++;
}
}
//计算每个像素在整幅图像中的比例
float maxPro = 0.0;
int kk = 0;
for (int 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 (int i = 0; i < GrayScale; i++)
{
w0 = w1 = u0tmp = u1tmp = u0 = u1 = u = deltaTmp = 0;
for (int 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;
u1 = u1tmp / w1;
u = u0tmp + u1tmp;
deltaTmp = w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);
if (deltaTmp > deltaMax)
{
deltaMax = deltaTmp;
threshold = i;
}
}
return threshold;
}
4.最大熵阈值分割算法
最大熵阈值分割(KSW)与Otsu算法类似,假设将图像分为背景和前景两个部分,熵代表信息量,图像信息量越大,熵就越大,最大熵法就是找出一个最佳阈值,来使得背景与前景两个部分熵之和最大。
实现:
//最大熵法阈值分割
int MaxEntropy_Segment(Mat&src,Mat&dst,int thresh,int p)
{
const int Grayscale = 256;
int Graynum[Grayscale]={0};
int height=src.rows;
int width=src.cols;
int i,j;
for(i=0;i<height;++i)
{
const uchar* ptr=src.ptr<uchar>(i);
for(j=0;j<width;++j)
{
if(ptr[j]==0) //排除掉图像中黑素像素点
continue;
Graynum[ptr[j]]++;
}
}
float probability=0.0;//概率
float max_Entropy;//最大熵
int totalpix=height*width;
for(i=0;i<Grayscale;++i)
{
float HO=0.0;//前景熵
float HB=0.0;//背景熵
//计算前景像素数
int frontpix=0;
for(j=0;j<i;++j)
{
frontpix+=Graynum[j];
}
//计算前景熵
for(int k=0;k<i;++k)
{
if(Graynum[k]!=0){
probability=(float)Graynum[k]/frontpix;
HO=HO+probability*log(1/probability);
}
}
//计算背景熵
for(int k=i;k<Grayscale;++k)
{
if(Graynum[k]!=0)
{
probability=(float)Graynum[k]/(totalpix-frontpix);
HB=HB+probability*log(1/probability);
}
}
//计算最大熵
if(HO+HB > max_Entropy)
{
max_Entropy=HO+HB;
thresh=i+p;
}
}
//阈值处理
src.copyTo(dst);
for(i=0;i<height;++i)
{
uchar* ptr=dst.ptr<ucahr>(i);
for(j=0;j<width;++j)
{
if(ptr[j]>thresh)
ptr[j]=255;
else
ptr[j]=0;
}
}
return thresh;
}