從貝葉斯理論到馬爾可夫隨機(jī)場(chǎng)--以圖像分割為例

• 馬爾可夫隨機(jī)場(chǎng)（CRF）
• 圖像分割過(guò)程
• Matlab代碼實(shí)現(xiàn)
• Python實(shí)現(xiàn)代碼
• 參考文獻(xiàn)

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

馬爾可夫隨機(jī)場(chǎng)（CRF）

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

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

圖像分割過(guò)程

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

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

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

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

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

先來(lái)說(shuō)說(shuō)P(W)這個(gè)標(biāo)記場(chǎng)w的先驗(yàn)概率。那么我們落實(shí)到圖像的每一個(gè)像素點(diǎn)上，也就是說(shuō)這個(gè)像素點(diǎn)是標(biāo)簽L的概率是多少，有人可能會(huì)說(shuō)，分割之前圖像的分類標(biāo)簽都不知道，還怎么算是某一類標(biāo)簽L的概率，這個(gè)問(wèn)題是有點(diǎn)較勁不好理解。

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

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

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

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

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

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

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

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

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

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

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

Matlab代碼實(shí)現(xiàn)

下面以條件迭代算法（ICM）來(lái)實(shí)例闡述MRF在圖像分割上的應(yīng)用，編程平臺(tái)為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對(duì)圖像進(jìn)行分割，類別數(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ù)分割');

• 基于馬爾可夫隨機(jī)場(chǎng)的圖像分割代碼

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

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

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

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

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

Python實(shí)現(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()

參考文獻(xiàn)

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

總結(jié)

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

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