本節主要介紹馬爾科夫的隨機場模型以及其用于圖像的分割算法中。基于馬爾科夫的隨機場（MRF）的圖像分割是一種基于統計的圖像分割算法，其模型參數少，空間約束性強，使用較為廣泛。

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

馬爾科夫隨機場在圖像領域的一大用途就是圖像的分割，本文也主要介紹該方法的圖像分割。

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

這是我們假設知道了最終分割結果W，但是W在開始是不知道的？現在的問題就是如何求W，我們所知道的不過是觀測到的圖像S，也就是說在知道S情況下求W，轉化為概率就是求取P(W|S) 問題就變成求取這個概率的最大值，也就是說根據S我們算出這個圖像的最有可能的分割標簽是什么。

基于此，問題才剛剛開始。我們由貝葉斯理論知道

觀察這個式子先不看能不能求出來，先給一些定義。因為W就是我們要求的分類標簽，我們稱P(W)為標記場w的先驗概率（要求它，事先又不知道，那就叫先驗概率吧，不知道是不是這樣來的）。P(S|W)是觀察值S的條件概率分布（也稱似然函數），看結構知道，這個是說在已知像素點標記w的情況下，它是真實的像素點s的概率，所以是一個似然函數，表示我的觀察像素點和真實的分類情況w像不像的一種程度。那P(S)是什么？S是觀察到的圖像，在分割前它就已經定了，我就要分割這幅圖像，那它還有概率嗎？沒有吧，所以P(S)認為是個定值。那么現在要求P(W|S)的最大值，也就是要求P(S|W)P(W)的最大值。

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

好，回到我們的問題，我們的問題要干嘛？要求一個像素點是標簽L的概率，開始標簽沒有，ok假設一個吧，一個不行，那么在開始我們把整幅圖像初始化一個分割標簽，盡管這個標簽不對（野雞），也可以假設一個稍微對一點的標簽（比如用其他分類算法（kmeans）得到一個預分割標簽）（這個時候認為是野雞100號）（好處在于這個時候算法不用迭代那么多次就可以停止了）。好了有了初始的標簽，那么再來求一個像素點是標簽L的概率。從這個時候就需要馬爾科夫協助了。我們想，要求一個像素點是標簽L的概率，此時這個像素點附近的領域都也已經劃分標簽了，雖然這個像素點也有一個標簽，但是這個標簽是上一代的標簽，那么它到下一代是需要進化的，馬爾科夫性質告訴我們這個像素點的分類情況只與附近一些領域分類情況有關，而與另一些領域沒有關系，也就是說我們可以根據這個像素點的附近領域的分類情況決定這個像素點在下一代是屬于哪一類的。

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

其中：

，是歸一化常數，參數 T 可控制P（w）的形狀，T越大越平坦。而，C為所有勢團集合，為勢團勢能，對

B為耦合系數，通常0.5-1。

糊不糊涂，不糊涂才怪，這么多關系，Gibbs也不容易，還沒完，那么什么是勢團，說白了就是和這個像素點構成一個小的檢測這個點連通性能的區域，一個圖如下：

總算完了，乍一看真是復雜，不知所云，其實細看還是可以理解的。說白了就是檢測一下這個點和周圍領域的相似程度。那些勢團就是要比較的對象。舉個例子，以一階勢團為例，假設某個點及其附近的8領域在某次迭代后的分類標簽是這樣的（假設分四類的情況）：

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

通觀Gibbs，其實就是一個求像素點與周圍像素點的不一致程度，或者說這個像素點的周圍像素點中，看看各個標簽出現的概率多大，就是這個點屬于對應標簽的概率。所以在實際編程上，也沒看到幾個人實實在在的用它的原版公式編程，反而簡化的計算，看看各個標簽出現的次數等等方式來計算的。

關于P(W)就說這么多，下面說說P(S|W)。P(S|W)，已知分類標簽那么它的像素值（灰度）是s的概率?，F在就假設w=1，某個像素點灰度為s，那么表示的意思就是在第一類里面像素灰度為s的概率。因為分類標簽在前面說到，每次迭代的時候是有一個分類標簽的（盡管不是最后的標簽），那么我們可以把屬于第一類的所有點都挑出來，考慮每個點都是獨立的，并且認為每一類里面的所有點服從高斯分布（正態分布），那么在每一類里面我們是不是可以根據這一類里面的這些點建立一個屬于這一類的高斯密度函數了，那么再來一個像素點值，把它帶到這類里面密度函數中去就可以得到這個概率了吧。同理對于2,3,4類，每一類都可以建立一個高斯密度函數，這樣就有四個高斯密度函數了，那么每一個點屬于這四類的概率就可以分別帶到這四個高斯密度函數中計算了。高斯密度函數一般形式為：

舉個例子假設現在我們求得的四類高斯函數，其圖像如下所示：

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

現在這一部分對于每個點也有了四個概率，結合上面每個點的四個概率，那么對每個點，屬于每一類的概率對應相乘，得到最終每個點屬于四個類的概率。這個時候就可以挑其中最大的那么就是該點所屬于的類了。這樣一次迭代，所有的點所屬于的類更新一遍，那這個新的類標簽作為下一次迭代的類標簽，反復下去，程序結束的條件可以是設置迭代次數，也可以是觀察類中心不在變化為止。

好了，上述的概率相乘取最大是原始一點的算法。其實可以看到，這個計算包括兩個部分，一般的講法，認為這兩部分每一部分都組成一個能量，換個說法就是最優化能量函數了，可以表示為如下：

像上述的概率相乘在實際中取一下對數就變成相加了。更深層次的馬爾科夫隨機場可以去看看相關論文，雖然方法可能不太一樣，但是原理上是一樣的。

下面以條件迭代算法（ICM）來實例闡述MRF在圖像分割上的應用，編程平臺為matlab，相關代碼如下：

clc

clear

close all

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

cluster_num = 4;%設置分類數

maxiter = 60;%最大迭代次數

%-------------隨機初始化標簽----------------

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

%-----------kmeans最為初始化預分割----------

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

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

iter = 0;while iter 《 maxiter

%-------計算先驗概率---------------

%這里我采用的是像素點和3*3領域的標簽相同

%與否來作為計算概率

%------收集上下左右斜等八個方向的標簽--------

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））;

%計算像素點8領域標簽相對于每一類的相同個數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;%計算概率

end

p_c（find（p_c == 0）） = 0.001;%防止出現0%---------------計算似然函數----------------

mu = zeros（1，4）;

sigma = zeros（1，4）;

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

index = find（label == i）;%找到每一類的點

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

%------計算每個像素點屬于每一類的似然概率--------

%------為了加速運算，將循環改為矩陣一起操作--------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

%找到聯合一起的最大概率最為標簽，取對數防止值太小

［~，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

將分類數設置為cluster_num=4，貼幾個中間結果：

上述程序除了注釋外再說一點，那就是對于是否需要預分類，程序里面寫了隨機分類和kmeans預分類兩者，貼的結果都是在隨機初始化的分類。當分別去實驗可以發現，如果采用kmeans預分類后，迭代開始分割的結果就很好，程序會在這個基礎上再變化一點，采用kmeans預分類的好處在于，分割的結果會更快達到穩定，且更準確。而采用隨機的預分類，分割的結果也算可以，這種方法得到的是一個局部最優解，當然kmeans得到的也是個局部最優解（不敢說最優解），但是前面的局部最優解比kmeans后的局部最優解又會差一點，尤其是在分割類別數大一點的時候，有kmeans預分類的效果會明顯好于隨機預分類。因為類別數越是大，相當于維度就越是大，那么它存在的局部最優解就越是多，那么從隨機的預分類（最差的情況）往最優解方向走時，中途可能隨便碰到一個局部最優解就走不動了。從某種意義上講，這個分割是一個np難問題，很難找到分割的最優解。