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

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

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

馬爾可夫隨機場（CRF）

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

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

圖像分割過程

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

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

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

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

W為要求的分類標簽，P(W)為標記場w的先驗概率，P(S|W)是觀察值S的條件概率分布，也稱為似然函數(shù)。在已知像素點標記w的情況下，真實的像素點的概率，所以是一個似然函數(shù)，所以P(S)認為是個定值。那么現(xiàn)在要求P(W|S)的最大值，也就是要求P(S|W)P(W)的最大值。

先來說說P(W)這個標記場w的先驗概率。那么我們落實到圖像的每一個像素點上，也就是說這個像素點是標簽L的概率是多少，有人可能會說，分割之前圖像的分類標簽都不知道，還怎么算是某一類標簽L的概率，這個問題是有點較勁不好理解。

我覺得，這就是一個雞蛋問題，是先有雞還是先有蛋的問題，我們要求這只雞，但是又沒有蛋怎么辦？一個簡單的辦法是我隨便弄一個蛋來，管他是雞蛋鴨蛋，有了蛋，雞（或者鴨）就可以有了，雖然這個雞不太像雞，比如說是只野雞，那么有了野雞，野雞再稍微進化一下，然后再下個蛋，蛋又長出野雞1號，野雞1號再稍微進化一下，然后再下個蛋，變成野雞2號，如此反復(fù)反復(fù)，知道野雞n號，然后驀然回首發(fā)現(xiàn)，野雞n號和我們所要的雞竟是那么的相似，那么把它就認為是我們要的雞了。有沒有糊涂？

其實優(yōu)化聚類里面很多問題都是雞蛋問題，比如曾經(jīng)介紹的FCM算法，有個先給u還是先給c的雞蛋問題，u的計算需要c，c的計算有需要u，那怎么辦，先假設(shè)一個吧。再比如EM算法的迭代過程也可以看成雞蛋問題，有了E步算M步，在回來算E步，在M步。。。。最終得到最優(yōu)解。

我們的問題要干嘛？要求一個像素點是標簽L的概率，開始標簽沒有，ok假設(shè)一個吧，一個不行，那么在開始我們把整幅圖像初始化一個分割標簽，盡管這個標簽不對（野雞），也可以假設(shè)一個稍微對一點的標簽（比如用其他分類算法（kmeans）得到一個預(yù)分割標簽）（這個時候認為是野雞100號）（好處在于這個時候算法不用迭代那么多次就可以停止了）。好了有了初始的標簽，那么再來求一個像素點是標簽L的概率。從這個時候就需要馬爾科夫協(xié)助了。

我們想，要求一個像素點是標簽L的概率，此時這個像素點附近的領(lǐng)域都也已經(jīng)劃分標簽了，雖然這個像素點也有一個標簽，但是這個標簽是上一代的標簽，那么它到下一代是需要進化的，馬爾科夫性質(zhì)告訴我們這個像素點的分類情況只與附近一些領(lǐng)域分類情況有關(guān)，而與另一些領(lǐng)域沒有關(guān)系，也就是說我們可以根據(jù)這個像素點的附近領(lǐng)域的分類情況決定這個像素點在下一代是屬于哪一類的。

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

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

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

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

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

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

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

Matlab代碼實現(xiàn)

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

• 首先將彩色圖像轉(zhuǎn)換為灰度圖像

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

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

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

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

clc clear close all img = double(imread('lena_gray.jpg')); %more quickly cluster_num = 4; %設(shè)置分類數(shù) maxiter = 60; %最大迭代次數(shù) %-------------隨機初始化標簽---------------- label = randi([1,cluster_num],size(img)); %-----------kmeans最為初始化預(yù)分割---------- % label = kmeans(img(:),cluster_num); % label = reshape(label,size(img)); iter = 0; while iter < maxiter%-------計算先驗概率---------------%這里我采用的是像素點和3*3領(lǐng)域的標簽相同%與否來作為計算概率%------收集上下左右斜等八個方向的標簽--------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領(lǐng)域標簽相對于每一類的相同個數(shù)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;%防止出現(xiàn)0%---------------計算似然函數(shù)----------------mu = zeros(1,4);sigma = zeros(1,4);%求出每一類的的高斯參數(shù)--均值方差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%------計算每個像素點屬于每一類的似然概率--------%------為了加速運算，將循環(huán)改為矩陣一起操作--------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 %找到聯(lián)合一起的最大概率最為標簽，取對數(shù)防止值太小[~,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

• 當(dāng)采用隨機初始化分割標簽，類別數(shù)為4類，迭代次數(shù)為60次時的分割結(jié)果如下。

• 當(dāng)采用K-Means初始化分割標簽，類別數(shù)為4類，迭代次數(shù)為60次時的分割結(jié)果如下。

• 當(dāng)采用K-Means初始化分割標簽，類別數(shù)為2類，迭代次數(shù)為60次時的分割結(jié)果如下。

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

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

Python實現(xiàn)代碼

• 將彩色圖片轉(zhuǎn)為灰度圖像

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。