圖像二值化的目的是最大限度的將圖象中感興趣的部分保留下來,在很多情況下,也是進行圖像分析、特征提取與模式識別之前的必要的圖像預處理過程。這個看似簡單的問題,在過去的四十年裡受到國內外學者的廣泛關注,產生了數以百計的閾值選取方法,但如同其他圖像分割算法一樣,沒有一個現有方法對各種各樣的圖像都能得到令人滿意的結果。
本文針對幾種經典而常用的二值發放進行了簡單的討論並給出了其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 '統計圖像的直方圖
Next
Next
For i = 0 To 255
Histgram(i) /= NumOfPixel '0到255的各灰度等級在圖像中各佔的比例
Next
For i = 0 To 255
Sum += Histgram(i)
If Sum >= Tile Then '得到閥值
Threshold = i
Exit For
End If
Next
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
Next
Next
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 '統計圖像的直方圖
Next
Next
For i = 0 To 255
AllSum += i * Histgram(i) ' 質量矩
AllPixelNumber += Histgram(i) ' 質量
Next
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
Next
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
Next
Next
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). 求出圖象的最大灰度值和最小灰度值,分別記為Zl和Zk,令初始閾值為:
(2). 根據閾值TK將圖象分割為前景和背景,分別求出兩者的平均灰度值Z0和ZB:
式中,Z(i,j)是圖像上(i,j)點的象素值,N(i,j)是(i,j)點的權值,一般取1。
(3)
(4). 若TK=TK+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)
Next
Next
NewThreshold = (MinGrayValue + MaxGrayValue) / 2
While Threshold <> NewThreshold And Iteration < 100
Threshold = NewThreshold
'根據閾值將圖像分割成目標和背景兩部分,求出兩部分的平均灰度值
For i = MinGrayValue To Threshold
IP1 += Histgram(i) * i
IS1 += Histgram(i)
Next
MeanGrayValue1 = CByte(IP1 / IS1)
For i = Threshold + 1 To MaxGrayValue
IP2 += Histgram(i) * i
IS2 += Histgram(i)
Next
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
Next
Next
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=H0 +HB,選取使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 '統計圖像的直方圖
Next
Next
For i = 0 To 255
Histgram(i) /= PixelNumber '統計各個灰度級出現的概率;
Next
'對每一個灰度級進行比較;
For i = 0 To 255
Temp1 = 0 : Temp2 = 0 : Temp = 0
For j = 0 To i
Temp += Histgram(j)
Next
For j = 0 To i
Temp1 += (-Histgram(j) / Temp * Log(Histgram(j) / Temp))
Next
For j = i + 1 To 255
Temp2 += (-Histgram(j) / (1 - Temp) * Log(Histgram(j) / (1 - Temp)))
Next
Sum(i) = Temp1 + Temp2
Next
MaxValue = 0.0 '找到使類的熵最大的灰度級;
For j = 0 To 255
If MaxValue < Sum(j) Then
MaxValue = Sum(j)
Threshold = j
End If
Next
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
Next
Next
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.mydown.com/news/2006-2-8/2785.html
我改寫成vb.net的了,可惜效果不理想,就沒有帖出效果圖了,不知道是不是我翻譯的不對。