參考：http://blog.csdn.net/on2way/article/details/47307927從貝葉斯理論到影象馬爾科夫隨機場

 劉偉強等；基於馬爾科夫隨機場的遙感影象分割和描述；東南大學學報；（29）：11-15，1999

version:2017.1.20 基本理解框架

Markov Random Filed

一、理論基礎

馬爾科夫隨機場是一種基於統計的影象分割演算法。它的特點主要為模型引數少、空間約束性強

馬爾科夫模型是一組事件的集合。在這個集合中，事件是逐個發生，並且下一刻事件的發生只由當前發生的事件決定，而與再之前的狀態沒有關係。拿博主【我愛i智慧】之天氣的例子說明為，【今天天氣好壞只與昨天天氣有關，而與前天乃至大前天都沒有關係】

二、基於影象分割應用

關於影象分割問題，從聚類角度講，就是一個影象聚類問題，把具有相同性質的畫素點設定為一類。也就是一個標籤分類問題，比如把一副影象分割成4類，那麼每一個畫素點必定屬於這四類中的某一類，假設四類為1，2，3，4類，那麼分割就是給每個畫素點找一個標籤類。好了，假設我們的待分割影象是S，大小m*n，那麼把每個畫素點放到集合S中，影象就是：，把待分割的影象稱為觀測到的影象。現在假設分割好了，每個畫素點都分了類，那麼把最終的分割結果稱為W，那麼顯然w大小與S一樣大，,其中所有的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的概率，此時這個畫素點附近的鄰域都已經劃分標籤了，雖然這個畫素點也有一個標籤，但是這個標籤是上一代的標籤，那麼它到下一代是需要進化的，馬爾科夫性質告訴我們這個畫素點的分類情況只與附近一些鄰域分類情況有關，而與另一些鄰域沒有關係，也就是說我們可以根據這個畫素點附近鄰域的分類情況決定這個畫素點在下一代是屬於哪一類的。

補充：馬爾科夫隨機場是以其區域性特性（馬爾科夫性）為特徵的，吉布斯隨機場是以其全域性特性（吉布斯分佈）為特徵的。Hammersley-Clifford定理則建立了這兩者之間的一致關係。

Hammersley-Clifford定理可知，要定義一個馬爾科夫隨機場，由於它與一個吉布斯隨機場相對應，如果定義了該吉布斯隨機場的能量函式，那麼這個馬爾科夫隨機場也就確定了

MAP：Maximum A Posterior，最大後驗概率估計

勢團勢能：（待補充）

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

其中，是歸一化常數，引數T可控制P（w）的形狀，T越大越平坦。而     ，C為所有勢團集合，  為勢團勢能。對  的定義如下： 

為耦合係數，通常0.5-1。 
糊不糊塗，不糊塗才怪，這麼多關係，Gibbs也不容易，還沒完，那麼什麼是勢團，說白了就是和這個畫素點構成一個小的【檢測這個點連通效能】的區域，一個圖如下： 

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

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

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

關於P(W)就說這麼多，下面說說P(S|W)。P(S|W)，已知分類標籤那麼它的畫素值（灰度）是s的概率。現在就假設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

1. 現在假設把影象分兩類cluster_num=2，採用隨機初始化分類，貼幾個之間結果： 
    
    
   

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

 


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