基于馬爾科夫的隨機(jī)場(chǎng)的圖像分割是一種基于統(tǒng)計(jì)的圖像分割算法

本節(jié)主要介紹馬爾科夫的隨機(jī)場(chǎng)模型以及其用于圖像的分割算法中?；隈R爾科夫的隨機(jī)場(chǎng)（MRF）的圖像分割是一種基于統(tǒng)計(jì)的圖像分割算法，其模型參數(shù)少，空間約束性強(qiáng)，使用較為廣泛。

首先了解一下馬爾科夫模型，純粹的馬爾科夫模型就是指一件事物的當(dāng)前狀態(tài)只與它之前的1個(gè)或者n個(gè)狀態(tài)有關(guān)，而與再之前的狀態(tài)沒(méi)有關(guān)系，比如今天天氣好壞只與昨天天氣有關(guān)，而與前天乃至大前天都沒(méi)有關(guān)系。符合這樣的一種特性的事物認(rèn)為其具有馬爾科夫性。那么引申到圖像領(lǐng)域，就是認(rèn)為圖像中某一點(diǎn)的特征（一般都是像素點(diǎn)灰色、顏色值等）只與其附近的一小塊領(lǐng)域有關(guān)，而與其他的領(lǐng)域無(wú)關(guān)。想想也是這樣，大多數(shù)時(shí)候，圖像中某一點(diǎn)像素和附近的像素是相關(guān)的，附近領(lǐng)域像素是黑的，那它八成也是黑的，很好理解。

馬爾科夫隨機(jī)場(chǎng)在圖像領(lǐng)域的一大用途就是圖像的分割，本文也主要介紹該方法的圖像分割。

關(guān)于圖像分割問(wèn)題，從聚類角度講，就是一個(gè)圖像的聚類問(wèn)題，把具有相同性質(zhì)的像素點(diǎn)設(shè)置為一類。在說(shuō)詳細(xì)點(diǎn)就是一個(gè)標(biāo)簽分類問(wèn)題，比如把一副圖像分割成4類，那么每一個(gè)像素點(diǎn)必定屬于這四類中的某一類，假設(shè)四類為1，2，3，4類，那么分割就是給每個(gè)像素點(diǎn)找一個(gè)標(biāo)簽類。好了，假設(shè)我們的待分割圖像是S，大小m*n，那么把每個(gè)像素點(diǎn)放到集合S中，圖像就是：S=S1,S2,...Sm*n把待分割的圖像稱為觀測(cè)到的圖像?，F(xiàn)在假設(shè)分割好了，每個(gè)像素點(diǎn)都分了類，那么把最終的分割結(jié)果稱為W，那么顯然w大小與S一樣大,W = W1,W2,...Wm*n, 其中所有的w取值都在1-L之間，L是最大的分類數(shù)，比如上面L=4，那么每個(gè)w都是在1-4之間取值的。

這是我們假設(shè)知道了最終分割結(jié)果W，但是W在開(kāi)始是不知道的？現(xiàn)在的問(wèn)題就是如何求W，我們所知道的不過(guò)是觀測(cè)到的圖像S，也就是說(shuō)在知道S情況下求W，轉(zhuǎn)化為概率就是求取P(W|S) 問(wèn)題就變成求取這個(gè)概率的最大值，也就是說(shuō)根據(jù)S我們算出這個(gè)圖像的最有可能的分割標(biāo)簽是什么。

基于此，問(wèn)題才剛剛開(kāi)始。我們由貝葉斯理論知道

觀察這個(gè)式子先不看能不能求出來(lái)，先給一些定義。因?yàn)閃就是我們要求的分類標(biāo)簽，我們稱P(W)為標(biāo)記場(chǎng)w的先驗(yàn)概率（要求它，事先又不知道，那就叫先驗(yàn)概率吧，不知道是不是這樣來(lái)的）。P(S|W)是觀察值S的條件概率分布（也稱似然函數(shù)），看結(jié)構(gòu)知道，這個(gè)是說(shuō)在已知像素點(diǎn)標(biāo)記w的情況下，它是真實(shí)的像素點(diǎn)s的概率，所以是一個(gè)似然函數(shù)，表示我的觀察像素點(diǎn)和真實(shí)的分類情況w像不像的一種程度。那P(S)是什么？S是觀察到的圖像，在分割前它就已經(jīng)定了，我就要分割這幅圖像，那它還有概率嗎？沒(méi)有吧，所以P(S)認(rèn)為是個(gè)定值。那么現(xiàn)在要求P(W|S)的最大值，也就是要求P(S|W)P(W)的最大值。

先來(lái)說(shuō)說(shuō)P(W)這個(gè)標(biāo)記場(chǎng)w的先驗(yàn)概率。那么我們落實(shí)到圖像的每一個(gè)像素點(diǎn)上，也就是說(shuō)這個(gè)像素點(diǎn)是標(biāo)簽L的概率是多少，有人可能會(huì)說(shuō)，分割之前圖像的分類標(biāo)簽都不知道，還怎么算是某一類標(biāo)簽L的概率，這個(gè)問(wèn)題是有點(diǎn)較勁不好理解。我覺(jué)得，這就是一個(gè)雞蛋問(wèn)題，是先有雞還是先有蛋的問(wèn)題，我們要求這只雞，但是又沒(méi)有蛋怎么辦？一個(gè)簡(jiǎn)單的辦法是我隨便弄一個(gè)蛋來(lái)，管他是雞蛋鴨蛋，有了蛋，雞（或者鴨）就可以有了，雖然這個(gè)雞不太像雞，比如說(shuō)是只野雞，那么有了野雞，野雞再稍微進(jìn)化一下，然后再下個(gè)蛋，蛋又長(zhǎng)出野雞1號(hào)，野雞1號(hào)再稍微進(jìn)化一下，然后再下個(gè)蛋，變成野雞2號(hào)，如此反復(fù)反復(fù)，知道野雞n號(hào)，然后驀然回首發(fā)現(xiàn)，野雞n號(hào)和我們所要的雞竟是那么的相似，那么把它就認(rèn)為是我們要的雞了。又沒(méi)有糊涂？其實(shí)優(yōu)化聚類里面很多問(wèn)題都是雞蛋問(wèn)題，比如曾經(jīng)介紹的FCM算法，有個(gè)先給u還是先給c的雞蛋問(wèn)題，u的計(jì)算需要c，c的計(jì)算有需要u，那怎么辦，先假設(shè)一個(gè)吧。再比如EM算法的迭代過(guò)程也可以看成雞蛋問(wèn)題，有了E步算M步，在回來(lái)算E步，在M步。。。。最終得到最優(yōu)解。

好，回到我們的問(wèn)題，我們的問(wèn)題要干嘛？要求一個(gè)像素點(diǎn)是標(biāo)簽L的概率，開(kāi)始標(biāo)簽沒(méi)有，ok假設(shè)一個(gè)吧，一個(gè)不行，那么在開(kāi)始我們把整幅圖像初始化一個(gè)分割標(biāo)簽，盡管這個(gè)標(biāo)簽不對(duì)（野雞），也可以假設(shè)一個(gè)稍微對(duì)一點(diǎn)的標(biāo)簽（比如用其他分類算法（kmeans）得到一個(gè)預(yù)分割標(biāo)簽）（這個(gè)時(shí)候認(rèn)為是野雞100號(hào)）（好處在于這個(gè)時(shí)候算法不用迭代那么多次就可以停止了）。好了有了初始的標(biāo)簽，那么再來(lái)求一個(gè)像素點(diǎn)是標(biāo)簽L的概率。從這個(gè)時(shí)候就需要馬爾科夫協(xié)助了。我們想，要求一個(gè)像素點(diǎn)是標(biāo)簽L的概率，此時(shí)這個(gè)像素點(diǎn)附近的領(lǐng)域都也已經(jīng)劃分標(biāo)簽了，雖然這個(gè)像素點(diǎn)也有一個(gè)標(biāo)簽，但是這個(gè)標(biāo)簽是上一代的標(biāo)簽，那么它到下一代是需要進(jìn)化的，馬爾科夫性質(zhì)告訴我們這個(gè)像素點(diǎn)的分類情況只與附近一些領(lǐng)域分類情況有關(guān)，而與另一些領(lǐng)域沒(méi)有關(guān)系，也就是說(shuō)我們可以根據(jù)這個(gè)像素點(diǎn)的附近領(lǐng)域的分類情況決定這個(gè)像素點(diǎn)在下一代是屬于哪一類的。

Ok既然我們認(rèn)為每個(gè)像素點(diǎn)的分類符合馬爾科夫隨機(jī)模型，而另一些學(xué)者證明了這種馬爾科夫隨機(jī)場(chǎng)可以與一個(gè)Gibbs隨機(jī)場(chǎng)等價(jià)（這就是Hammcrslcy-Clifford定理，人家證明的，那就這樣叫吧）。而Gibbs隨機(jī)場(chǎng)也有一個(gè)概率密度函數(shù)，這樣我們就用求圖像的Gibbs隨機(jī)場(chǎng)的概率P代替了P(W)。那么Gibbs隨機(jī)場(chǎng)的概率P怎么表示呢？如下：

其中：

，是歸一化常數(shù)，參數(shù) T 可控制P（w）的形狀，T越大越平坦。而，C為所有勢(shì)團(tuán)集合，為勢(shì)團(tuán)勢(shì)能，對(duì)

B為耦合系數(shù)，通常0.5-1。

糊不糊涂，不糊涂才怪，這么多關(guān)系，Gibbs也不容易，還沒(méi)完，那么什么是勢(shì)團(tuán)，說(shuō)白了就是和這個(gè)像素點(diǎn)構(gòu)成一個(gè)小的檢測(cè)這個(gè)點(diǎn)連通性能的區(qū)域，一個(gè)圖如下：

總算完了，乍一看真是復(fù)雜，不知所云，其實(shí)細(xì)看還是可以理解的。說(shuō)白了就是檢測(cè)一下這個(gè)點(diǎn)和周圍領(lǐng)域的相似程度。那些勢(shì)團(tuán)就是要比較的對(duì)象。舉個(gè)例子，以一階勢(shì)團(tuán)為例，假設(shè)某個(gè)點(diǎn)及其附近的8領(lǐng)域在某次迭代后的分類標(biāo)簽是這樣的（假設(shè)分四類的情況）：

現(xiàn)在要計(jì)算中心的像素點(diǎn)在下一次迭代中是屬于第幾類（這一代是第3類），ok采用一階勢(shì)能，這里需要說(shuō)明一點(diǎn)，這個(gè)像素點(diǎn)無(wú)非是1-4之間的一類，那么我們需要分別計(jì)算下一代它是第1,2,3,4類時(shí)的勢(shì)團(tuán)。先假設(shè)這個(gè)點(diǎn)是第1類，先比較左右，發(fā)現(xiàn)都是1-1一樣，ok記一下B,在與上B,與下，不一樣，那么記一下-B，如果再來(lái)勢(shì)團(tuán)的話，斜上方，1-2，不一樣-B，斜下方，1-1一樣B,一次類推，就可以將中心點(diǎn)如果是第1類的勢(shì)團(tuán)計(jì)算出來(lái)。那么在計(jì)算中心點(diǎn)如果是第2類，發(fā)現(xiàn)這個(gè)時(shí)候除了3個(gè)斜著的方向是一樣的外，其他的都不一樣。在計(jì)算假設(shè)是第三類，第四類。這樣有了每類下的勢(shì)團(tuán)，就可以計(jì)算概率了，根據(jù)公式往回帶呀帶，最終得到這個(gè)點(diǎn)如果屬于第一類的概率，第二類的概率，第三類的概率，第四類的概率，這四個(gè)值在后面會(huì)一起用到來(lái)決定這個(gè)點(diǎn)到底屬于哪一類。你可能會(huì)說(shuō)，知道了這四個(gè)概率，看哪個(gè)最大不就出來(lái)了嗎？注意這還只是第一項(xiàng)的先驗(yàn)概率，后面還有個(gè)似然函數(shù)的概率沒(méi)有計(jì)算呢，后面你會(huì)看到這個(gè)似然函數(shù)的概率對(duì)于每個(gè)像素點(diǎn)也有四個(gè)概率，在這四個(gè)概率都計(jì)算出來(lái)之后，那么這個(gè)點(diǎn)屬于第一類的概率用著兩個(gè)都屬于第一類的概率相乘，第二類也是如此，最后再比較這四個(gè)概率的最大值作為要標(biāo)記的標(biāo)簽。

通觀Gibbs，其實(shí)就是一個(gè)求像素點(diǎn)與周圍像素點(diǎn)的不一致程度，或者說(shuō)這個(gè)像素點(diǎn)的周圍像素點(diǎn)中，看看各個(gè)標(biāo)簽出現(xiàn)的概率多大，就是這個(gè)點(diǎn)屬于對(duì)應(yīng)標(biāo)簽的概率。所以在實(shí)際編程上，也沒(méi)看到幾個(gè)人實(shí)實(shí)在在的用它的原版公式編程，反而簡(jiǎn)化的計(jì)算，看看各個(gè)標(biāo)簽出現(xiàn)的次數(shù)等等方式來(lái)計(jì)算的。

關(guān)于P(W)就說(shuō)這么多，下面說(shuō)說(shuō)P(S|W)。P(S|W)，已知分類標(biāo)簽?zāi)敲此南袼刂担ɑ叶龋┦莝的概率?，F(xiàn)在就假設(shè)w=1，某個(gè)像素點(diǎn)灰度為s，那么表示的意思就是在第一類里面像素灰度為s的概率。因?yàn)榉诸悩?biāo)簽在前面說(shuō)到，每次迭代的時(shí)候是有一個(gè)分類標(biāo)簽的（盡管不是最后的標(biāo)簽），那么我們可以把屬于第一類的所有點(diǎn)都挑出來(lái)，考慮每個(gè)點(diǎn)都是獨(dú)立的，并且認(rèn)為每一類里面的所有點(diǎn)服從高斯分布（正態(tài)分布），那么在每一類里面我們是不是可以根據(jù)這一類里面的這些點(diǎn)建立一個(gè)屬于這一類的高斯密度函數(shù)了，那么再來(lái)一個(gè)像素點(diǎn)值，把它帶到這類里面密度函數(shù)中去就可以得到這個(gè)概率了吧。同理對(duì)于2,3,4類，每一類都可以建立一個(gè)高斯密度函數(shù)，這樣就有四個(gè)高斯密度函數(shù)了，那么每一個(gè)點(diǎn)屬于這四類的概率就可以分別帶到這四個(gè)高斯密度函數(shù)中計(jì)算了。高斯密度函數(shù)一般形式為：

舉個(gè)例子假設(shè)現(xiàn)在我們求得的四類高斯函數(shù)，其圖像如下所示：

某一點(diǎn)的灰度為s=110，那么它對(duì)應(yīng)的分別P(s|W1)、P(s|W2)、P(s|W3)、P(s|W3)，就分別如上圖所示呢，很顯然的可以看到110在第三類下的概率最大，其他的幾個(gè)概率也可以出來(lái)的。

現(xiàn)在這一部分對(duì)于每個(gè)點(diǎn)也有了四個(gè)概率，結(jié)合上面每個(gè)點(diǎn)的四個(gè)概率，那么對(duì)每個(gè)點(diǎn)，屬于每一類的概率對(duì)應(yīng)相乘，得到最終每個(gè)點(diǎn)屬于四個(gè)類的概率。這個(gè)時(shí)候就可以挑其中最大的那么就是該點(diǎn)所屬于的類了。這樣一次迭代，所有的點(diǎn)所屬于的類更新一遍，那這個(gè)新的類標(biāo)簽作為下一次迭代的類標(biāo)簽，反復(fù)下去，程序結(jié)束的條件可以是設(shè)置迭代次數(shù)，也可以是觀察類中心不在變化為止。

好了，上述的概率相乘取最大是原始一點(diǎn)的算法。其實(shí)可以看到，這個(gè)計(jì)算包括兩個(gè)部分，一般的講法，認(rèn)為這兩部分每一部分都組成一個(gè)能量，換個(gè)說(shuō)法就是最優(yōu)化能量函數(shù)了，可以表示為如下：

像上述的概率相乘在實(shí)際中取一下對(duì)數(shù)就變成相加了。更深層次的馬爾科夫隨機(jī)場(chǎng)可以去看看相關(guān)論文，雖然方法可能不太一樣，但是原理上是一樣的。

下面以條件迭代算法（ICM）來(lái)實(shí)例闡述MRF在圖像分割上的應(yīng)用，編程平臺(tái)為matlab，相關(guān)代碼如下：

clc

clear

close all

img = double（imread（‘lena.jpg’））;%more quickly

cluster_num = 4;%設(shè)置分類數(shù)

maxiter = 60;%最大迭代次數(shù)

%-------------隨機(jī)初始化標(biāo)簽----------------

label = randi（［1，cluster_num］，size（img））;

%-----------kmeans最為初始化預(yù)分割----------

% label = kmeans（img（：），cluster_num）;

% label = reshape（label，size（img））;

iter = 0;while iter 《 maxiter

%-------計(jì)算先驗(yàn)概率---------------

%這里我采用的是像素點(diǎn)和3*3領(lǐng)域的標(biāo)簽相同

%與否來(lái)作為計(jì)算概率

%------收集上下左右斜等八個(gè)方向的標(biāo)簽--------

label_u = imfilter（label，［0，1，0;0，0，0;0，0，0］，‘replicate’）;

label_d = imfilter（label，［0，0，0;0，0，0;0，1，0］，‘replicate’）;

label_l = imfilter（label，［0，0，0;1，0，0;0，0，0］，‘replicate’）;

label_r = imfilter（label，［0，0，0;0，0，1;0，0，0］，‘replicate’）;

label_ul = imfilter（label，［1，0，0;0，0，0;0，0，0］，‘replicate’）;

label_ur = imfilter（label，［0，0，1;0，0，0;0，0，0］，‘replicate’）;

label_dl = imfilter（label，［0，0，0;0，0，0;1，0，0］，‘replicate’）;

label_dr = imfilter（label，［0，0，0;0，0，0;0，0，1］，‘replicate’）;

p_c = zeros（4，size（label，1）*size（label，2））;

%計(jì)算像素點(diǎn)8領(lǐng)域標(biāo)簽相對(duì)于每一類的相同個(gè)數(shù)for i = 1:cluster_num

label_i = i * ones（size（label））;

temp = ~（label_i - label_u） + ~（label_i - label_d） + 。。.

~（label_i - label_l） + ~（label_i - label_r） + 。。.

~（label_i - label_ul） + ~（label_i - label_ur） + 。。.

~（label_i - label_dl） +~（label_i - label_dr）;

p_c（i，：） = temp（：）/8;%計(jì)算概率

end

p_c（find（p_c == 0）） = 0.001;%防止出現(xiàn)0%---------------計(jì)算似然函數(shù)----------------

mu = zeros（1，4）;

sigma = zeros（1，4）;

%求出每一類的的高斯參數(shù)--均值方差for i = 1:cluster_num

index = find（label == i）;%找到每一類的點(diǎn)

data_c = img（index）;

mu（i） = mean（data_c）;%均值

sigma（i） = var（data_c）;%方差

end

p_sc = zeros（4，size（label，1）*size（label，2））;

% for i = 1:size（img，1）*size（img，2）

% for j = 1:cluster_num

% p_sc（j，i） = 1/sqrt（2*pi*sigma（j））*.。。

% exp（-（img（i）-mu（j））^2/2/sigma（j））;

% end

% end

%------計(jì)算每個(gè)像素點(diǎn)屬于每一類的似然概率--------

%------為了加速運(yùn)算，將循環(huán)改為矩陣一起操作--------for j = 1:cluster_num

MU = repmat（mu（j），size（img，1）*size（img，2），1）;

p_sc（j，：） = 1/sqrt（2*pi*sigma（j））*.。。

exp（-（img（：）-MU）.^2/2/sigma（j））;

end

%找到聯(lián)合一起的最大概率最為標(biāo)簽，取對(duì)數(shù)防止值太小

［~，label］ = max（log（p_c） + log（p_sc））;

%改大小便于顯示

label = reshape（label，size（img））;

%---------顯示----------------if ~mod（iter，6）

figure;

n=1;

end

subplot（2，3，n）;

imshow（label，［］）

title（［‘iter = ’，num2str（iter）］）;

pause（0.1）;

n = n+1;

iter = iter + 1;

end

將分類數(shù)設(shè)置為cluster_num=4，貼幾個(gè)中間結(jié)果：

上述程序除了注釋外再說(shuō)一點(diǎn)，那就是對(duì)于是否需要預(yù)分類，程序里面寫了隨機(jī)分類和kmeans預(yù)分類兩者，貼的結(jié)果都是在隨機(jī)初始化的分類。當(dāng)分別去實(shí)驗(yàn)可以發(fā)現(xiàn)，如果采用kmeans預(yù)分類后，迭代開(kāi)始分割的結(jié)果就很好，程序會(huì)在這個(gè)基礎(chǔ)上再變化一點(diǎn)，采用kmeans預(yù)分類的好處在于，分割的結(jié)果會(huì)更快達(dá)到穩(wěn)定，且更準(zhǔn)確。而采用隨機(jī)的預(yù)分類，分割的結(jié)果也算可以，這種方法得到的是一個(gè)局部最優(yōu)解，當(dāng)然kmeans得到的也是個(gè)局部最優(yōu)解（不敢說(shuō)最優(yōu)解），但是前面的局部最優(yōu)解比kmeans后的局部最優(yōu)解又會(huì)差一點(diǎn)，尤其是在分割類別數(shù)大一點(diǎn)的時(shí)候，有kmeans預(yù)分類的效果會(huì)明顯好于隨機(jī)預(yù)分類。因?yàn)轭悇e數(shù)越是大，相當(dāng)于維度就越是大，那么它存在的局部最優(yōu)解就越是多，那么從隨機(jī)的預(yù)分類（最差的情況）往最優(yōu)解方向走時(shí)，中途可能隨便碰到一個(gè)局部最優(yōu)解就走不動(dòng)了。從某種意義上講，這個(gè)分割是一個(gè)np難問(wèn)題，很難找到分割的最優(yōu)解。