從貝葉斯理論到馬爾可夫隨機場--以圖像分割為例

• 馬爾可夫隨機場（CRF）
• 圖像分割過程
• Matlab代碼實現
• Python實現代碼
• 參考文獻

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

馬爾可夫隨機場（CRF）

首先了解以下馬爾可夫模型，純粹的馬爾可夫模型是指一個事物的當前狀態只與之前的1個狀態有關，而與其他狀態沒有關系。比如今天的天氣好壞只與昨天的天氣有關，而與前天以及大前天都沒有關系，符合這樣的一種特性的事物認為其具有馬爾可夫性。

那么引申到圖像領域，就是認為圖像中某一點的特征（一般都是像素點的灰度值或顏色值）只與其附近的一小塊鄰域有關，而與其他區域無關。

圖像分割過程

馬爾可夫隨機場在圖像領域的重要應用就是圖像分割。關于圖像分割問題，從聚類角度講就是一個圖像聚類問題，把具有相同性質的像素點聚集為一類。簡而言之，圖像分類就是像素級別的分類問題，比如把一個圖像分割為四個部分，就是把所有的像素點分為四類，即給每個像素點找一個標簽類。

假設待分割的圖像是S，大小為m*n，那么把每個像素點放到集合S中，圖像就是$S=\{S_1,S_2,?,S_{m?n}\}$，把待分割的圖像稱為觀測圖像。假設圖像已經分割完成，每個像素都有具體的類別，把最終的分割結果稱為W，顯然W的shape與S的shape相同，每個類別的取值范圍是1-4.

現在的問題就是如何求W，我們所知道就是觀測到的圖像S，也就是說知道S的情況下求W，轉化為概率就是求取$P(W|S)$，圖像分割問題就變成了求取這個概率的最大值，也就是說根據S計算出圖像的分割標簽。

根據貝葉斯公式$$P(Y|X)=\frac{P(X|Y)P(Y)}{P(X)}$$.

W為要求的分類標簽，P(W)為標記場w的先驗概率，P(S|W)是觀察值S的條件概率分布，也稱為似然函數。在已知像素點標記w的情況下，真實的像素點的概率，所以是一個似然函數，所以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的概率，此時這個像素點附近的領域都也已經劃分標簽了，雖然這個像素點也有一個標簽，但是這個標簽是上一代的標簽，那么它到下一代是需要進化的，馬爾科夫性質告訴我們這個像素點的分類情況只與附近一些領域分類情況有關，而與另一些領域沒有關系，也就是說我們可以根據這個像素點的附近領域的分類情況決定這個像素點在下一代是屬于哪一類的。

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

$$P(w)=\frac{1}{z}exp(?\frac{1}{T}{U}_2(w))$$

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

$$P(x|w_i)=\frac{1}{\sqrt{2\pi }\sigma }exp(?\frac{(x?u{)}^2}{2{\sigma }^2})$$

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

某一點的灰度為s=110，那么它對應的分別P(s|W1)、P(s|W2)、P(s|W3)、P(s|W4)，就分別如上圖所示呢，很顯然的可以看到110在第二類下的概率最大，其他的幾個概率也可以出來的。
現在這一部分對于每個點也有了四個概率，結合上面每個點的四個概率，那么對每個點，屬于每一類的概率對應相乘，得到最終每個點屬于四個類的概率。這個時候就可以挑其中最大的那么就是該點所屬于的類了。

這樣一次迭代，所有的點所屬于的類更新一遍，那這個新的類標簽作為下一次迭代的類標簽，反復下去，程序結束的條件可以是設置迭代次數，也可以是觀察類中心不在變化為止。
好了，上述的概率相乘取最大是原始一點的算法。其實可以看到，這個計算包括兩個部分，一般的講法，認為這兩部分每一部分都組成一個能量，換個說法就是最優化能量函數了，可以表示為如下：
$$W=argmin(U_1(w,S),U_2(w))$$
像上述的概率相乘在實際中取一下對數就變成相加了。更深層次的馬爾科夫隨機場可以去看看相關論文，雖然方法可能不太一樣，但是原理上是一樣的。

Matlab代碼實現

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

• 首先將彩色圖像轉換為灰度圖像

img_rgb=imread('lena.jpg'); %讀入需要轉換的彩色圖像 imshow(img_rgb); %顯示彩色圖像 img_gray=rgb2gray(img_rgb); %RGB圖像轉為灰度圖像 figure imshow(img_gray); %顯示灰度圖像 imwrite(img_gray,'lena_gray.jpg'); %保存灰度圖像

• 使用K-Means對圖像進行分割，類別數為2

clc clear close all img = double(imread('lena_gray.jpg')); %more quickly cluster_num = 2; %設置分類數 %-----------kmeans初始化預分割---------- label = kmeans(img(:),cluster_num); label = reshape(label,size(img)); imshow(label,[]); title('kmeans初始化預分割');

• 基于馬爾可夫隨機場的圖像分割代碼

clc clear close all img = double(imread('lena_gray.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_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(find(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;iter = iter + 1; end

• 當采用隨機初始化分割標簽，類別數為4類，迭代次數為60次時的分割結果如下。

• 當采用K-Means初始化分割標簽，類別數為4類，迭代次數為60次時的分割結果如下。

• 當采用K-Means初始化分割標簽，類別數為2類，迭代次數為60次時的分割結果如下。

上述程序除了注釋外再說一點，那就是對于是否需要預分類，程序里面寫了隨機分類和kmeans預分類兩者。當分別去實驗可以發現，如果采用kmeans預分類后，迭代開始分割的結果就很好，程序會在這個基礎上再變化一點。采用kmeans預分類的好處在于，分割的結果會更快達到穩定，且更準確。

而采用隨機的預分類，分割的結果也算可以，這種方法得到的是一個局部最優解，當然kmeans得到的也是個局部最優解（不敢說最優解），但是前面的局部最優解比kmeans后的局部最優解又會差一點，尤其是在分割類別數大一點的時候，有kmeans預分類的效果會明顯好于隨機預分類。因為類別數越是大，相當于維度就越是大，那么它存在的局部最優解就越是多，那么從隨機的預分類（最差的情況）往最優解方向走時，中途可能隨便碰到一個局部最優解就走不動了。從某種意義上講，這個分割是一個np難問題，很難找到分割的最優解。

Python實現代碼

• 將彩色圖片轉為灰度圖像

import cv2 import matplotlib.pyplot as pltimg = cv2.imread('./dataset/images/lena.jpg', cv2.IMREAD_COLOR) img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB) print('RGB:', img.shape) # RGB: (512, 512, 3) plt.imshow(img) plt.xticks([]) plt.yticks([]) plt.show()img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY) print('GRAY:', img.shape) # GRAY: (512, 512) plt.imshow(img, cmap='gray') plt.xticks([]) plt.yticks([]) plt.show()

參考文獻

https://blog.csdn.net/on2way/article/details/47307927

總結

以上是生活随笔為你收集整理的从贝叶斯理论到马尔可夫随机场（MRF）--以图像分割为例的全部內容，希望文章能夠幫你解決所遇到的問題。

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