几种经典的二值化方法及其vb.net实现

梦想的边缘 2008-05-05

展开全文

几种经典的二值化方法及其vb.net实现

http://blog.csdn.net/laviewpbt/archive/2006/05/24/753774.aspx

图像二值化的目的是最大限度的将图象中感兴趣的部分保留下来，在很多情况下，也是进行图像分析、特征提取与模式识别之前的必要的图像预处理过程。这个看似简单的问题，在过去的四十年里受到国内外学者的广泛关注，产生了数以百计的阈值选取方法，但如同其他图像分割算法一样，没有一个现有方法对各种各样的图像都能得到令人满意的结果。

本文针对几种经典而常用的二值发放进行了简单的讨论并给出了其vb.net 实现。

1、P-Tile法

Doyle于1962年提出的P-Tile (即P分位数法)可以说是最古老的一种阈值选取方法。该方法根据先验概率来设定阈值，使得二值化后的目标或背景像素比例等于先验概率，该方法简单高效，但是对于先验概率难于估计的图像却无能为力。
这样，我们先统计出图像中各灰度等级象素出现的总次数，接着得到对应的概率，找到一个合适的阈值，使得小于该阈值的象素数大于或等于先验概率。具体实现如下：

'程序实现功能：经典的二值化方法及其实现
'作者： laviewpbt
'联系方式： laviewpbt@sina.com
'QQ:33184777
'版本：Version 1.0.0
'说明：复制请保留源作者信息，转载请说明，欢迎大家提出意见和建议

' ******************************************************************************************

' 函数名： P_Tile

' 功能： P分位数法二值化图像

' 参数： Bmp ------ Bitmap 待处理位图

' Tile ------ Single 先验概率

' TimeElapse ------ Integer 处理所需的时间

' 返回值： Boolean

' 作者： laviewpbt

' 时间： 2005-5-20 12：48

' 修改者：

' 修改时间：

' ******************************************************************************************

Public Shared Function P_Tile(ByVal Bmp As Bitmap, ByVal Tile As Single, ByVal TimeElapse As Integer) As Byte

If Tile < 0 OrElse Tile > 1 Then Throw New Exception("Tile的值只可以在0和1之间")

TimeElapse = Environment.TickCount

Dim i, j, Stride, Temp As Integer

Dim NumOfPixel As Integer = Bmp.Width * Bmp.Height ‘总象素数

Dim Histgram(255), Sum As Single

Dim BmpData(), Threshold As Byte

ReadBitmap(Bmp, BmpData) '读取数据

Stride = (((Bmp.Width * 24) + 31) \ 32) * 4

For i = 0 To Bmp.Height - 1

For j = 0 To Bmp.Width - 1

Histgram(BmpData(i * Stride + j * 3)) += 1 '统计图像的直方图

For i = 0 To 255

Histgram(i) /= NumOfPixel '0到255的各灰度等级在图像中各占的比例

For i = 0 To 255

Sum += Histgram(i)

If Sum >= Tile Then '得到阀值

Threshold = i

Exit For

End If

For i = 0 To Bmp.Height - 1

For j = 0 To Bmp.Width - 1

Temp = i * Stride + j * 3 '二值显示

If BmpData(Temp) <= Threshold Then

BmpData(Temp) = 0 : BmpData(Temp + 1) = 0 : BmpData(Temp + 2) = 0

Else

BmpData(Temp) = 255 : BmpData(Temp + 1) = 255 : BmpData(Temp + 2) = 255

End If

WriteBitmap(Bmp, BmpData) '写入数据

TimeElapse = Environment.TickCount - TimeElapse

Return Threshold

End Function

如果在上述程序中设置Tile为0.5，则整个二值化的图片中有黑白各占一半左右。

原始图像

Tile 取0.8时的结果，用时0豪秒,阈值180

在原图中，地面绿色的区域差不多占了整个图片的80%，而Tile 取0.8时的分割效果也还比较理想。

2、OTSU 算法（大津法）

OSTU算法可以说是自适应计算单阈值（用来转换灰度图像为二值图像）的简单高效方法。1978 OTSU年提出的最大类间方差法以其计算简单、稳定有效，一直广为使用。

' ******************************************************************************************

' 函数名： Otsu

' 功能： Otsu二值图像

' 参数： Bmp ------ Bitmap 待处理位图

' TimeElapse ------ Integer 处理所需的时间

' 返回值： Boolean

' 作者： laviewpbt

' 时间： 2005-5-20 3：45

' 修改者：

' 修改时间：

' ******************************************************************************************

Public Shared Function Otsu(ByVal Bmp As Bitmap, ByRef TimeElapse As Integer) As Byte

TimeElapse = Environment.TickCount

Dim i, j, k, Stride, Temp As Integer

Dim AllSum, SumSmall, SumBig, PartSum As Integer

Dim AllPixelNumber, PixelNumberSmall, PixelNumberBig As Integer

Dim ProbabilitySmall, ProbabilityBig, Probability, MaxValue As Double

Dim BmpData(), Threshold As Byte

Dim Histgram(255) As Integer '图像直方图，256个点

Dim Width As Integer = Bmp.Width, Height As Integer = Bmp.Height

Dim PixelNumber As Integer = Bmp.Width * Bmp.Height

Stride = (((Bmp.Width * 24) + 31) \ 32) * 4

ReadBitmap(Bmp, BmpData)

Dim Number As Integer = BmpData.Length - 1

For i = 0 To Bmp.Height - 1

For j = 0 To Bmp.Width - 1

Histgram(BmpData(i * Stride + j * 3)) += 1 '统计图像的直方图

For i = 0 To 255

AllSum += i * Histgram(i) ' 质量矩

AllPixelNumber += Histgram(i) ' 质量

MaxValue = -1.0

For i = 0 To 255

PixelNumberSmall += Histgram(i)

PixelNumberBig = AllPixelNumber - PixelNumberSmall

If PixelNumberBig = 0 Then Exit For

SumSmall += i * Histgram(i)

SumBig = AllSum - SumSmall

ProbabilitySmall = CDbl(SumSmall) / PixelNumberSmall

ProbabilityBig = CDbl(SumBig) / PixelNumberBig

Probability = PixelNumberSmall * ProbabilitySmall * ProbabilitySmall + PixelNumberBig * ProbabilityBig * ProbabilityBig

If Probability > MaxValue Then

MaxValue = Probability

Threshold = i

End If

For i = 0 To Bmp.Height - 1

For j = 0 To Bmp.Width - 1

Temp = i * Stride + j * 3 '二值显示

If BmpData(Temp) <= Threshold Then

BmpData(Temp) = 0 : BmpData(Temp + 1) = 0 : BmpData(Temp + 2) = 0

Else

BmpData(Temp) = 255 : BmpData(Temp + 1) = 255 : BmpData(Temp + 2) = 255

End If

WriteBitmap(Bmp, BmpData)

TimeElapse = Environment.TickCount - TimeElapse

Return Threshold

End Function

把上述过程中的 Probability = PixelNumberSmall * ProbabilitySmall * ProbabilitySmall + PixelNumberBig * ProbabilityBig * ProbabilityBig 改为 Probability = PixelNumberSmall * PixelNumberBig * (ProbabilityBig - ProbabilitySmall) * (ProbabilityBig - ProbabilitySmall) 改为也能得到较为合理的结果。

大津法选取出来的阈值非常理想，对各种情况的表现都较为良好。虽然它在很多情况下都不是最佳的分割，但分割质量通常都有一定的保障，可以说是最稳定的分割。

Probability 取第一中参数的时的效果（上部白色和背景重叠了），用时50ms,阈值123

Probability 取第二种参数的时的效果,用时52ms,阈值248

3、迭代法（最佳阀值法）

(1)．求出图象的最大灰度值和最小灰度值，分别记为Z_l和Z_k，令初始阈值为：

(2)．根据阈值TK将图象分割为前景和背景，分别求出两者的平均灰度值Z₀和Z_B:

式中，Z(i,j）是图像上(i,j)点的象素值，N(i,j)是(i,j)点的权值，一般取1。

(3)

    (4)．  若T^K=T^K+1，则所得即为阈值，否则转2，迭代计算。

    以下给出迭代求阈值的实现：

' ******************************************************************************************

' 函数名： BestThreshold

' 功能：最佳阀值法二值图像

' 参数： Bmp ------ Bitmap 待处理位图

' TimeElapse ------ Integer 处理所需的时间

' 返回值： Boolean

' 作者： laviewpbt

' 时间： 2005-5-20 4：20

' 修改者：

' 修改时间：

' ******************************************************************************************

Public Shared Function BestThreshold(ByVal Bmp As Bitmap, ByRef TimeElapse As Integer) As Byte

TimeElapse = Environment.TickCount

Dim i, j, k, Stride, Temp As Integer

Dim Threshold ,NewThreshold,MaxGrayValue ,MinGrayValue,MeanGrayValue1,MeanGrayValue2 As Byte

Dim IP1 As Integer, IP2 As Integer, IS1 As Integer, IS2 As Integer

Dim Iteration As Integer, Histgram(255) As Integer

MaxGrayValue = 0 : MinGrayValue = 255

Stride = (((Bmp.Width * 24) + 31) \ 32) * 4

Dim BmpData() As Byte

ReadBitmap(Bmp, BmpData)

Dim Number As Integer = BmpData.Length - 1

'求出图像中的最小和最大灰度值，并计算阈值初值为

For i = 0 To Bmp.Height - 1

For j = 0 To Bmp.Width - 1

Temp = i * Stride + j * 3

Histgram(BmpData(Temp)) += 1 '统计图像的直方图

If MinGrayValue > BmpData(Temp) Then MinGrayValue = BmpData(Temp)

If MaxGrayValue < BmpData(Temp) Then MaxGrayValue = BmpData(Temp)

NewThreshold = (MinGrayValue + MaxGrayValue) / 2

While Threshold <> NewThreshold And Iteration < 100

Threshold = NewThreshold

'根据阈值将图像分割成目标和背景两部分，求出两部分的平均灰度值

For i = MinGrayValue To Threshold

IP1 += Histgram(i) * i

IS1 += Histgram(i)

MeanGrayValue1 = CByte(IP1 / IS1)

For i = Threshold + 1 To MaxGrayValue

IP2 += Histgram(i) * i

IS2 += Histgram(i)

MeanGrayValue2 = CByte(IP2 / IS2)

'求出新的阈值：

NewThreshold = (MinGrayValue + MaxGrayValue) / 2

Iteration += 1

End While

For i = 0 To Bmp.Height - 1

For j = 0 To Bmp.Width - 1

Temp = i * Stride + j * 3 '二值显示

If BmpData(Temp) <= Threshold Then

BmpData(Temp) = 0 : BmpData(Temp + 1) = 0 : BmpData(Temp + 2) = 0

Else

BmpData(Temp) = 255 : BmpData(Temp + 1) = 255 : BmpData(Temp + 2) = 255

End If

WriteBitmap(Bmp, BmpData)

TimeElapse = Environment.TickCount - TimeElapse

Return True

End Function

最佳阀值法效果，用时60ms,阈值128

4、一维最大熵阈值法

它的思想是统计图像中每一个灰度级出现的概率，计算该灰度级的熵 ,假设以灰度级T分割图像，图像中低于T灰度级的像素点构成目标物体（O），高于灰度级T的像素点构成背景（B），那么各个灰度级在本区的分布概率为：

　　O区： i=1,2……,t

　　B区： i=t+1,t+2……L-1

　　上式中的 ,这样对于数字图像中的目标和背景区域的熵分别为：
　　

对图像中的每一个灰度级分别求取W=H₀ +H_B，选取使W最大的灰度级作为分割图像的阈值，这就是一维最大熵阈值图像分割法。

' ******************************************************************************************

' 函数名： MaxEntropy

' 功能：一维最大熵二值化图像

' 参数： Bmp ------ Bitmap 待处理位图

' TimeElapse ------ Integer 处理所需的时间

' 返回值： Boolean

' 作者： laviewpbt

' 时间： 2005-5-20 5：40

' 修改者：

' 修改时间：

' ******************************************************************************************

Public Shared Function MaxEntropy(ByVal Bmp As Bitmap, ByRef TimeElapse As Integer) As Byte

TimeElapse = Environment.TickCount

Dim i, j, Stride, Threshold As Integer

Dim Temp1, Temp2, Temp, Sum(255), MaxValue As Single

Dim BmpData() As Byte

Dim Histgram(255) As Single

Dim Width As Integer = Bmp.Width, Height As Integer = Bmp.Height

Stride = (((Bmp.Width * 24) + 31) \ 32) * 4

Dim PixelNumber As Integer = Bmp.Width * Bmp.Height

ReadBitmap(Bmp, BmpData)

Dim Number As Integer = BmpData.Length - 1

For i = 0 To Bmp.Height - 1

For j = 0 To Bmp.Width - 1

Histgram(BmpData(i * Stride + j * 3)) += 1 '统计图像的直方图

For i = 0 To 255

Histgram(i) /= PixelNumber '统计各个灰度级出现的概率；

'对每一个灰度级进行比较；

For i = 0 To 255

Temp1 = 0 : Temp2 = 0 : Temp = 0

For j = 0 To i

Temp += Histgram(j)

For j = 0 To i

Temp1 += (-Histgram(j) / Temp * Log(Histgram(j) / Temp))

For j = i + 1 To 255

Temp2 += (-Histgram(j) / (1 - Temp) * Log(Histgram(j) / (1 - Temp)))

Sum(i) = Temp1 + Temp2

MaxValue = 0.0 '找到使类的熵最大的灰度级；

For j = 0 To 255

If MaxValue < Sum(j) Then

MaxValue = Sum(j)

Threshold = j

End If

For i = 0 To Bmp.Height - 1

For j = 0 To Bmp.Width - 1

Temp = i * Stride + j * 3 '二值显示

If BmpData(Temp) <= Threshold Then

BmpData(Temp) = 0 : BmpData(Temp + 1) = 0 : BmpData(Temp + 2) = 0

Else

BmpData(Temp) = 255 : BmpData(Temp + 1) = 255 : BmpData(Temp + 2) = 255

End If

WriteBitmap(Bmp, BmpData)

TimeElapse = Environment.TickCount - TimeElapse

Return Threshold

End Function

一维最大熵效果，用时150ms,阈值124

这种方法的缺点是仅仅考虑了像素点的灰度信息，没有考虑到像素点的空间信息，所以当图像的信噪比降低时分割效果不理想。不过二维最大熵法我还没有搞定，呵呵。

5 、聚类算法

聚类算法是把一副图像分割成n个类（n>=2），当n=2时也可以作为二值化的一种有效算法，其基本思想是把某一象素点归纳入距离其最近的一类中，通过不段迭代，直到两次迭带的结果符合指定的精度为止。一般步骤如下：

(1)指定类别数n、初始聚类中心，迭代停止参数theta;

(2)计算每点到到各类的距离，并将改点归入距离其最近的一类，计算新的聚类中心；

(3) 比较两次聚类中心的差异程度是否小于theat，是则停止迭带，否转第二步。

上述过程中距离的概念是广义的，除了我们常用的欧式距离外，还可以取其他的能够描述两类差异的参数，如相关系数。但要注意不同参数距离最近的意义不同，比如如果取欧式距离，则数值越小，距离越近，而取相关系数时，数值越大，距离越近。

除了上述硬聚类的方法，把模糊理论运用到聚类中的方法也得到了广泛的应用，因为模糊现象更加符合自然界的规律，著名的模糊聚类算法有FCM，模糊K均值等，但模糊算法比硬聚类的方法所需要的计算时间长，在实际中可以把硬聚类的中心作为模糊聚类的初始中心，一提高运算的实时性。

对于上述图片，运用HCM和FCM算法的二值化效果如下（实际上HCM或FCM分割后并不是黑白图片，而是具有两种彩色的图片，这里为了比较把其中一种颜色显示为黑色，一种为白色）。

HCM二值化的效果

FCM二值化的效果

在初始中心方面，我取的是[0，0，0]，[255，255，255]，HCM用时大约0.6s左右，FCM用时大约2s(这里没有采用HCM的最终聚类中心，也是[0，0，0]，[255，255，255]）。所以相对来说速度慢了不少，不过嘛，HCM本来不是主要干这个的，但看得出聚类的方法得到的结果比较细致。

了解了聚类的过程，改写成代码不是很复杂的过程（呵呵，现在看来不复杂，我最开始写的时候可是用了一个星期，那时真菜）。详细代码见：

模糊聚类算法（FCM）和硬聚类算法（HCM)的VB6.0实现及其应用

6、区域生长法

有关于该算法的代码请参看http://www./news/2006-2-8/2785.html

我改写成vb.net的了，可惜效果不理想，就没有帖出效果图了，不知道是不是我翻译的不对。