图像二值化阈值自动选取


原来找的些东西,先放在这里,有时间再整理,
http://www.programfan.com/article/showarticle.asp?id=2219 图像分割是图像处理这门学科中的基础难题,基于阈值的分割则又是图像分割的最基本的难题之一,其难点在于阈值的选取。事实证明,阈值的选择的恰当与否对分割的效果起着决定性的作用。由于阈值选取对图像分割的基础性,本文主要在【1】、【2】、【3】、【4】等的基础上,对一些当前流行的阈值选取算法做了探讨、实现和比较。多阈值分割虽然能进一步提高图像分割的质量,但由于它只是分割技巧的处理问题,而与单阈值分割并无本质的区别。因此本文并不对多阈值分割进行讨论,而只考虑单阈值分割的情形。

1. 双峰法

双峰法的原理及其简单:它认为图像由前景和背景组成,在灰度直方图上,前后二景都形成高峰,在双峰之间的最低谷处就是图像的阈值所在。根据这一原理,我们给出了它的实现,部分代码如下(Pascal语言描述,以下同):

//intPeak、intPeak2、intValley:峰值和直方图值 


//intIndx::相应的灰度值 


intPeak,intIndx,intPeak2,intIndx2,intValley,intValleyIndx:integer; 


//初始双峰值 


intPeak:=0; 


intPeak2:=0; 




//取得第一峰值 


for intLoop:=0 to 255 do 


if intPeak<=intGrayLevel[intLoop] then 


begin 


intPeak:=intGrayLevel[intLoop]; 


intIndx:=intLoop; 


end; 




//取得第二峰值 


for intLoop:=0 to 255 do 


Begin 


if (intPeak2<=intGrayLevel[intLoop]) and (intLoop<>intIndx) then 


begin 


intPeak2:=intGrayLevel[intLoop]; 


intIndx2:=intLoop; 


end 


end; 




//取得双峰之间的谷值 


intValley:=intSize; 


if intIndx2 
 
 

for intLoop:=intIndx2 to intIndx do 
 

if intValley>intGrayLevel[intLoop] then 
 

begin 
 

intValley:=intGrayLevel[intLoop]; 
 

intValleyIndx:=intLoop; 
 

end;



从分割的效果来看,当前后景的对比较为强烈时,分割效果较好;否则基本无效。

2. 迭代法

迭代法是基于逼近的思想,其步骤如下:

1. 求出图象的最大灰度值和最小灰度值,分别记为ZMAX和ZMIN,令初始阈值T0=(ZMAX+ZMIN)/2;

2. 根据阈值TK将图象分割为前景和背景,分别求出两者的平均灰度值ZO和ZB;

3. 求出新阈值TK+1=(ZO+ZB)/2;

4. 若TK=TK+1,则所得即为阈值;否则转2,迭代计算。

以下给出迭代求阈值的部分实现:

//阈值初始为0 
 

intThresholdVal:=0; 
 

intThresholdVal2:=0; 
 



//总灰度值 
 

intTotalGrayLevel:=0; 
 

for intLoop:=0 to 255 do 
 

if intGrayLevel[intLoop]<>0 then 
 

intTotalGrayLevel:=intTotalGrayLevel+intLoop*intGrayLevel[intLoop]; 
 



//求出初始最大灰度值 
 

for intLoop:=0 to 255 do 
 

if intGrayLevel[intLoop]>0 then 
 

begin 
 

intLGrayLevel:=intLoop; 
 

intThresholdVal:=intLoop; 
 

break; 
 

end; 
 



//求出初始最小灰度值和初始阈值 
 

for intLoop:=255 downto 0 do 
 

if intGrayLevel[intLoop]>0 then 
 

begin 
 

intRGrayLevel:=intLoop; 
 

intThresholdVal:=(intThresholdVal+intLoop)div 2; 
 

break; 
 

end; 
 



//迭代求解 
 

while intThresholdVal<>intThresholdVal2 do 
 

begin 
 

intThresholdVal2:=intThresholdVal; 
 

intCount:=0; 
 

intLGrayLevel:=0; 
 

for intLoop:=0 to intThresholdVal do 
 

if intGrayLevel[intLoop]<>0 then 
 

begin 
 

intCount:=intCount+intGrayLevel[intLoop]; 
 

intLGrayLevel:=intLGrayLevel+intLoop*intGrayLevel[intLoop]; 
 

end; 
 

intRGrayLevel:=intTotalGrayLevel-intLGrayLevel; 
 

intLGrayLevel:=intLGrayLevel div intCount; 
 

intRGrayLevel:=intRGrayLevel div (intSize-intCount); 
 

intThresholdVal:=(intLGrayLevel+intRGrayLevel)div 2; 
 

end;



迭代所得的阈值分割的图象效果良好。基于迭代的阈值能区分出图像的前景和背景的主要区域所在,但在图像的细微处(如图1中的浅色线条)还没有很好的区分度。

但令人惊讶的是,对某些特定图象,微小数据的变化却会引起分割效果的巨大改变,两者的数据只是稍有变化,但分割效果却反差极大,个中原因还有待进一步研究。

3. 大津法(OTSU法)

大津法由大津于1979年提出,对图像Image,记t为前景与背景的分割阈值,前景点数占图像比例为w0, 平均灰度为u0;背景点数占图像比例为w1,平均灰度为u1。图像的总平均灰度为:u=w0*u0+w1*u1。从最小灰度值到最大灰度值遍历t,当t使得值g=w0*(u0-u)2+w1*(u1-u)2 最大时t即为分割的最佳阈值。对大津法可作如下理解:该式实际上就是类间方差值,阈值t分割出的前景和背景两部分构成了整幅图像,而前景取值u0,概率为w0,背景取值u1,概率为w1,总均值为u,根据方差的定义即得该式。因方差是灰度分布均匀性的一种度量,方差值越大,说明构成图像的两部分差别越大,当部分目标错分为背景或部分背景错分为目标都会导致两部分差别变小,因此使类间方差最大的分割意味着错分概率最小。

直接应用大津法计算量较大,因此我们在实现时采用了等价的公式g=w0*w1*(u0-u1)2。部分计算过程如下:


//遍历所有灰度值求Max g。 
 

for intCurrentLevel:=0 to intArrLen do 
 

begin 
 

if intSclGrayLevel[intCurrentLevel]=0 then 
 

continue 
 

else 
 

begin 
 

//计算当阈值为intCurrentLevel时的g 
 

intCount:=0; 
 

intSumPels:=0; 
 

for intLoop:=0 to intCurrentLevel do 
 

begin 
 

intCount:=intCount+intSclGrayLevel[intLoop]; 
 

intSumPels:=intSumPels+intSumPelsArr[intLoop]; 
 

end; 
 

w0:=intCount/intSize; 
 

u0:=intSumPels/intCount; 
 

w1:=1-w0; 
 

if intSize-intCount<>0 then 
 

u1:=(intTotalPels-intSumPels)/(intSize-intCount) 
 

else 
 

u1:=0; 
 



RlTempO:=w0*w1*(u0-u1)*(u0-u1); 
 

if RlTempO>RlMaxO then 
 

begin 
 

RlMaxO:=RlTempO; 
 

Result:=intCurrentLevel; 
 

end; 
 

end;



我们在测试中发现:大津法选取出来的阈值非常理想,对各种情况的表现都较为良好。虽然它在很多情况下都不是最佳的分割,但分割质量通常都有一定的保障,可以说是最稳定的分割。由上可知,大津算法是一种较为通用的分割算法。在它的思想的启迪下,人们进一步提出了多种类似的评估阈值的算法,具体可参加【5】、【6】等。

4. 灰度拉伸-一种改进的大津法

大津法得到了广泛的应用,但有人发现,大津法致命的缺陷是当目标物与背景灰度差不明显时,会出现无法忍受的大块黑色区域,甚至会丢失整幅图像的信息。为了解决这个问题,有人提出了灰度拉伸的增强大津法。这种方法的原理其实就是在大津法的基础上通过增加灰度的级数来增强前后景的灰度差,从而解决问题。灰度增加的方法是用原有的灰度级乘上同一个系数,从而扩大灰度的级数,特别地,当乘上的系数为1时,这就是大津法的原型,因此,大津法可以看做是这种方法的一个特例。

在实现中,我们实现了多种灰度拉伸,发现对不同的图像,当遇上不同的拉伸系数时,分割效果也相差甚远。

5. Kirsh算子

在【4】中提出了基于Kirsh算子的分割方法,其思想为:对数字图像的每个像素i,考虑它的八个邻点的灰度值,以其中三个相邻点的加权和减去剩下五个邻点的加权和得到差值,令三个邻点绕该像素点不断移位,取此八个差值的最大值作为Kirsh算子。即:设Si为三邻点之和,Ti为五邻点之和,则Kirsh算子定义为K(i)=max{1,max〔5Si-3Ti〕}如取阈值THk,则当K(i)>THk时,像素i为阶跃边缘点。此外,【4】的作者认为:假设图像大小为H×W个像素点,其边缘点像素一般不会超过5×H个。基于这一假设,该文作者提出:(对一幅图像)用Kirsh算法,取某一较低的初始阈值THk(以保证目标和背景间灰度变化很小的图像边缘也能被取出),对于每个像素点i计算其Kirsh算子,如果K(i)>THk,则i为边缘点,边缘点数N(初始值为0)加1,一旦边缘点数超过5×H ,而i还小于整幅图像的像素数,说明阈值取得太低,致使许多不是边缘点的像素也被取出,因此需提高阈值。如此反复,即可获得分割图像所需的阈值。

但在实现中,本文作者发现,【4】中的叙述颇有值得探讨之处,如在H×W图像中,H和W之间的关系是完全对称的,两者之间如何抉择?此外,在求Kirsh算子K(i)=max{1,max〔5Si-3Ti〕}时也颇有疑虑之处,由其求得的结果分割图像效果并不明显。基于对称性和归一化的考虑,笔者把Kirsh算子改为:K(i)=max{1,max abs(5Si-3Ti) div 15 },并根据在实际运行中的效果,对W和H的选取为:if W>H then use 5*H else use 5*W。在实际应用中表明,修改后的分割质量显著提高。但与【4】文中作者声称的效果及其示例相比,仍有相当的距离,特别是它不能解决前后景对比不强烈时的分割情形。但当前后背景对比十分强烈且集中时,Kirsh算子法却会有十分突出的表现。
http://forum.assuredigit.com/display_topic_threads.asp?ForumID=8&TopicID=3480
/*
OTSU 算法可以说是自适应计算单阈值(用来转换灰度图像为二值图像)的简单高效方法。下面的代码最早由 Ryan Dibble提供,此后经过多人Joerg.Schulenburg, R.Z.Liu 等修改,补正。

算法对输入的灰度图像的直方图进行分析,将直方图分成两个部分,使得两部分之间的距离最大。划分点就是求得的阈值。

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 (unsigned char *image, int rows, int cols, int x0, int y0, int dx, int dy, int vvv) 
 
{ 
 
unsigned char *np; // 图像指针 
 
int thresholdValue=1; // 阈值 
 
int ihist[256]; // 图像直方图,256个点 
 
int i, j, k; // various counters 
 
int n, n1, n2, gmin, gmax; 
 
double m1, m2, sum, csum, fmax, sb; 
 
// 对直方图置零... 
 
memset(ihist, 0, sizeof(ihist)); 
 
gmin=255; gmax=0; 
 
// 生成直方图 
 
for (i = y0 + 1; i < y0 + dy - 1; i++) { 
 
 np = &image[i*cols+x0+1]; 
 
for (j = x0 + 1; j < x0 + dx - 1; j++) { 
 
ihist[*np]++; 
 
if(*np > gmax) gmax=*np; 
 
if(*np < gmin) gmin=*np; 
 
 np++; /* next pixel */ 
 
} 
 
} 
 
// set up everything 
 
sum = csum = 0.0; 
 
n = 0; 
 
for (k = 0; k <= 255; k++) { 
 
sum += (double) k * (double) ihist[k]; /* x*f(x) 质量矩*/ 
 
n += ihist[k]; /* f(x) 质量 */ 
 
} 
 
if (!n) { 
 
// if n has no value, there is problems... 
 
fprintf (stderr, "NOT NORMAL thresholdValue = 160 "); 
 
return (160); 
 
} 
 
// do the otsu global thresholding method 
 
fmax = -1.0; 
 
n1 = 0; 
 
for (k = 0; k < 255; k++) { 
 
n1 += ihist[k]; 
 
if (!n1) { continue; } 
 
n2 = n - n1; 
 
if (n2 == 0) { break; } 
 
csum += (double) k *ihist[k]; 
 
m1 = csum / n1; 
 
m2 = (sum - csum) / n2; 
 
sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2); 
 
/* bbg: note: can be optimized. */ 
 
if (sb > fmax) { 
 
fmax = sb; 
 
thresholdValue = k; 
 
} 
 
} 
 
// at this point we have our thresholding value 
 
// debug code to display thresholding values 
 
if ( vvv & 1 ) 
 
fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d gmax=%d ", 
 
thresholdValue, gmin, gmax); 
 
return(thresholdValue); 
 
}