由于經常有讀者說運行出錯，我又本地校準了下程序，由于版權限制，lena圖MATLAB新版已經撤除了，這里改成了cameraman的圖

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

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

version:2017.1.20 基本理解框架

Markov Random Filed

一、理論基礎

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

馬爾科夫模型是一組事件的集合。在這個集合中，事件是逐個發生，并且下一刻事件的發生只由當前發生的事件決定，而與再之前的狀態沒有關系。拿博主【我愛i智能】之天氣的例子說明為，【今天天氣好壞只與昨天天氣有關，而與前天乃至大前天都沒有關系】，引申到圖像領域為：圖像中某一點特征（一般指灰度值、RGB值等）只與其附近的一小塊鄰域有關，而與其他的鄰域無關。

二、基于圖像分割應用

關于圖像分割問題，從聚類角度講，就是一個圖像聚類問題，把具有相同性質的像素點設置為一類。也就是一個標簽分類問題，比如把一副圖像分割成4類，那么每一個像素點必定屬于這四類中的某一類，假設四類為1，2，3，4類，那么分割就是給每個像素點找一個標簽類。好了，假設我們的待分割圖像是S，大小m*n，那么把每個像素點放到集合S中，圖像就是：，把待分割的圖像稱為觀測到的圖像?，F在假設分割好了，每個像素點都分了類，那么把最終的分割結果稱為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怎么表示呢？如下：

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

?

?

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

通觀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('cameraman.tif'));%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)); for iter=0:maxiter-1%-------計算先驗概率---------------%這里采用的是像素點和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_numlabel_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;%計算概率endp_c((p_c == 0)) = 0.001;%防止出現0%---------------計算似然函數----------------mu = zeros(1,4);sigma = zeros(1,4);%求出每一類的的高斯參數--均值方差for i = 1:cluster_numindex = find(label == i);%找到每一類的點data_c = img(index);mu(i) = mean(data_c);%均值sigma(i) = var(data_c);%方差endp_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_numMU = 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;endsubplot(2,3,n);imshow(label,[])title(['iter = ',num2str(iter)]);pause(0.1);n = n+1; end

現在假設把圖像分兩類cluster_num=2，采用隨機初始化分類，貼幾個之間結果：?
?
?

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

?


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

總結

以上是生活随笔為你收集整理的马尔科夫随机场之图像分割【二】的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。