MRF馬爾可夫隨機場入門
Intro
MRF是一種普遍應用於圖像分割的模型，固然我看到MRF的時候並非由於分割，而是在圖像生成領域，有的paper利用MRF模型來生成圖像，所以入門一下MRF，並以分割模型爲例記一下代碼。python
Model
Target
在圖像分割中，咱們的任務是給定一張圖像，輸出每一個像素的標籤。所以咱們就是要獲得在給定圖片特徵之下，標籤機率最大化時所對應的標籤。dom
所以能夠這麼建模： $$ \hat{\omega} = arg \max_{\omega \in \Omega} P(\omega|f) $$ 其中w表示標籤，f表示圖像特徵，求最大後驗機率。函數
根據貝葉斯理論，上式右邊能夠寫成： $$ P(\omega|f) = \frac{P(f|\omega)P(\omega)}{P(f)} $$ 其中，P(f)是常量，由於當一張圖片肯定以後，P(f)便肯定了。所以，上式只取決於分子部分。分子又能夠表達爲$P(f,\omega)$，因此咱們直接建模的實際上是這個部分，計算的也是這個部分，這是與CRF不一樣的一點(MRF是直接對左邊建模，不分解爲右邊，因此沒個樣本都要算一遍後驗機率，而後乘起來最大化，MRF實際上是經過對等式右邊分子建模"曲線救國")。優化
所以，咱們的任務中只須要對分子的兩個部分進行定義便可。ui
Neighbors
像素Neighbors的定義很簡單，就是這個像素周圍的其餘像素。spa
舉例而言，下圖分別是中心點像素的四鄰域和八鄰域。翻譯
Hammersley-Clifford Theorem
定理的內容爲：code
若是一個分佈$P(x)>0$知足無向圖$G$中的局部馬爾可夫性質，當且僅當$P(x)$能夠表示爲一系列定義在最大團上的非負函數的乘積形式，即： $$ P(x) = \frac{1}{Z} \prod_{c \in C} \phi (x_c) $$ 其中$C$爲$G$中最大團集合，也就是全部的最大團組成的集合，$\phi(x_c) \ge 0$是定義在團$c$上的勢能函數，Z是配分函數，用來將乘積歸一化爲機率的形式。 $$ Z = \sum_{x \in \Chi } \prod_{c \in C} \phi (x_c) $$ 無向圖模型與有向圖模型的一個重要區別就是配分函數Z。orm
Hammersley-Clifford Theorem代表，無向圖模型和吉布斯分佈是一致的，因此將$P(\omega)$定義下式： $$ P(\omega) = \frac{1}{Z}exp(-U(\omega)) = \frac{1}{Z}exp(- \sum_{c \in C} V_c(\omega)) $$ 其中，Z做爲normalization項，$Z = \sum exp(-U(\omega))$，U定義爲勢能，而等號最右邊將U變成了V的求和，在後面咱們會說到，這裏實際上是每一個原子團的勢能的求和。blog
Clique
Clique就是咱們上面提到的「團」的概念。集合$c$是$S$的原子團當且僅當c中的每一個元素都與該集合中的其餘元素相鄰。那麼Clique就是全部$c$的並集。 $$ C = c_1 \cup c_2 \cdots c_n $$
舉例而言：
一個像素的四鄰域及他本身組成的集合的原子團能夠分爲singleton和doubleton如圖所示。
Clique Potential
翻譯過來就是勢能，用$V(w)$表示，描述的是一個Clique的能量。
那麼，一個像素的領域的勢能就是每一個團的能量的和。 $$ U(\omega) = \sum_{c \in C} V_c(w) $$ 其中c表示原子團,c表示Clique，V是如何定義的呢？
在圖像分割中，能夠以一階團爲例， $$ V_c(\omega) = \beta \delta(w,w_s) = \left{\begin{aligned}\beta &&w = w_s \-\beta && w \neq w_s \\end{aligned}\right. $$ 到這裏，$P(\omega)$的全部變量解釋完了，下一步是計算$P(f|\omega)$
$P(f|\omega)$的計算
$P(f|\omega)$被認爲是服從高斯分佈的，也就是說，若是咱們知道了這個像素的標籤是什麼，那麼他的像素值應該服從這個標籤下的條件機率的高斯分佈。其實他服從高斯分佈仍是很好理解的，咱們已知這個像素點的label好比說是A，那麼咱們去統計一下全部標籤是A的點的像素值的均值和方差，顯然以這個均值和方差爲參數的高斯分佈更加契合這裏的條件分佈。 $$ P(f_s|\omega_s) = \frac{1}{\sqrt{2\pi}\sigma_{w_s}}exp(-\frac{(f_s - \mu_{\omega_s})^2}{2\sigma^2_{\omega_s}}) $$ 計算每一個類別的像素均值和方差，帶入公式，即得條件機率。
最後，就是最大化$P(\omega)P(f|\omega)$,以對數形式轉化爲求和的形式去優化，最大化$log(P(\omega)) + log(P(f|\omega))$.
Coding
import numpy as np
import cv2 as cv
import copy
class MRF():
def __init__(self,img,max_iter = 100,num_clusters = 5,init_func = None,beta = 10):
self.max_iter = max_iter
self.kernels = np.zeros(shape = (8,3,3))
self.beta = beta
self.num_clusters = num_clusters
for i in range(9):
if i < 4:
self.kernels[i,i//3,i%3] = 1
elif i > 4:
self.kernels[i-1,(i-1)//3,(i-1)%3] = 1
self.img = img
if init_func is None:
self.labels = np.random.randint(low = 1,high = num_clusters + 1,size = img.shape,dtype = np.uint8)
def __call__(self):
#print(self.labels)
#print(self.img)
img = self.img.reshape((-1,))
for iter in range(self.max_iter):
p1 = np.zeros(shape = (self.num_clusters,self.img.shape[0] * self.img.shape[1]))
for cluster_idx in range(self.num_clusters):
temp = np.zeros(shape = (self.img.shape))
for i in range(8):
res = cv.filter2D(self.labels,-1,self.kernels[i,:,:])
temp[(res == (cluster_idx + 1))] += self.beta
temp[(res == (cluster_idx + 1))] -= self.beta
#temp += (res == (cluster_idx + 1)).astype(np.int)
temp = np.exp(-temp)
#temp = temp/8.0
p1[cluster_idx,:] = temp.reshape((-1,))
p1 = p1 / np.sum(p1,axis = 0)
p1[p1 == 0] = 1e-3
mu = np.zeros(shape = (self.num_clusters,))
sigma = np.zeros(shape = (self.num_clusters,))
for i in range(self.num_clusters):
#mu[i] = np.mean(self.img[self.labels == (i+1)])
data = self.img[self.labels == (i+1)]
if np.sum(data) > 0:
mu[i] = np.mean(data)
sigma[i] = np.var(data)
else:
mu[i]= 0
sigma[i] = 1
#print(sigma[i])
#sigma[sigma == 0] = 1e-3
p2 = np.zeros(shape = (self.num_clusters,self.img.shape[0] * self.img.shape[1]))
for i in range(self.img.shape[0] * self.img.shape[1]):
for j in range(self.num_clusters):
#print(sigma[j])
p2[j,i] = -np.log(np.sqrt(2*np.pi)*sigma[j]) -(img[i]-mu[j])**2/2/sigma[j];
self.labels = np.argmax(np.log(p1) + p2,axis = 0) + 1
self.labels = np.reshape(self.labels,self.img.shape).astype(np.uint8)
print("-----------start-----------")
print(p1)
print("-" * 20)
print(p2)
print("----------end------------")
#print("iter {} over!".format(iter))
#self.show()
#print(self.labels)
def show(self):
img = copy.copy(self.labels)
img[img == 1] = 0
img[img == 2] = 80
img[img == 3] = 160
img[img == 4] = 240
#img = self.labels / (self.num_clusters) * 255
cv.imshow("res",img.astype(np.uint8))
cv.waitKey(0)
if __name__ == "__main__":
img = cv.imread("/home/xueaoru/圖片/test002.jpg")
img = cv.cvtColor(img,cv.COLOR_BGR2GRAY)
img = img/255.
#img = np.random.rand(64,64)
#img = cv.resize(img,(256,256))
mrf = MRF(img = img,max_iter = 20,num_clusters = 2)
mrf()
mrf.show()
#print(mrf.kernels)
Input:
Output(num_clusters = 4):
Output(num_clusters = 2):