机器学习24：概率图--马尔可夫随机场(MRF)
机器学习24：概率图--马尔可夫随机场(MRF)
1，马尔可夫随机场(MarkovRandom Field)：
下图是一个简单的马尔可夫随机场：
马尔可夫随机场是典型的马尔可夫网，图中每个结点表示一个或一组变量，结点之间的边表示两个变量之间的依赖关系。
马尔可夫随机场是一种具有马尔可夫性的随机场，要理解什么是马尔可夫随机场，我们得要先理解什么是随机场，而什么又是马尔可夫性。
1.1，随机场：
在概率论中, 由样本空间Ω = {0, 1, …, G − 1}n取样构成的随机变量Xi所组成的S = {X1, …, Xn}。若对所有的ω∈Ω，式子π(ω) > 0均成立，则称π为一个随机场。
当给每一个位置中按照某种分布随机赋予相空间的一个值之后，其全体就叫做随机场。我们不妨拿种地来打个比方。其中有两个概念：位置（site），相空间（phase space）。“位置”好比是一亩亩农田；“相空间”好比是种的各种庄稼。我们可以给不同的地种上不同的庄稼，这就好比给随机场的每个“位置”，赋予相空间里不同的值。所以，俗气点说，随机场就是在哪块地里种什么庄稼的事情。
1.2，团，极大团：
团（Clique）：图中节点的子集，其中任意两个节点之间都有边连接；极大团：一个团，其中加入任何一个其他的节点都不能再形成团。
1.3，HammersleyClifford 定理：
马尔可夫随机场中，多个变量之间的联合概率分布可以基于团分解为多个势函数的乘积，每个势函数仅与一个团相关。
比如，对于下图所表达的马尔科夫随机场：
它的表达就是：
1.4，分离集：
设 A、B、C 都是马尔可夫随机场中的节点集合，若从 A 中的节点到 B 中的节点都必须经过 C 中的 节点，则称 A 和 B 被 C 分离，C 称为分离集（Separating Set）。如下图：
为了便于讨论，将上图中的A, B, C分别对应于单变量xA,xB,xC，于是上图可以简化为：
于是，对于上图所代表的马尔可夫随机场，由HammersleyClifford 定理可得联合概率：
1.5，条件独立性：
下面证明，在给定xC时，xA与xB条件独立：
因此：
1.6，马尔可夫性：
全局马尔可夫性：设结点集合A，B是在无向图G中被结点集合C分开的任意结点集合，则在给定随机变量xC的条件下，随机变量xA和xB条件独立。
由全局马尔可夫性可以得到两个有用的推论：
局部马尔可夫性：设v是无向图G中任意一个结点，W是与v有边相连的所有结点，G中其他结点记做O；则在给定随机变量xW的条件下，随机变量xV和xO条件独立。此时这里的W是v的分离集，局部马尔可夫性是全局马尔可夫性的推论。
成对马尔可夫性：设u和v是无向图G中任意两个没有边直接连接的结点，G中其他结点的集合记做O；则在给定随机变量xO的条件下，随机变量xU和xV条件独立。
满足上述马尔可夫性，三者其一的无向图称之为MRF。
1.7，势函数的选择：
马尔可夫随机场有一组势函数(potentialfunctions)，也称为因子。势函数是定义在变量子集上的非负实函数，主要用于定义概率分布函数。
显然，势函数ΨQ(XQ)的作用是定量刻画变量集XQ（一个极大团）中变量之间的相关关系，它应该是非负函数，且在所偏好的变量取值上有较大的函数值。比如：假定XQ中的变量均是二值变量，若势函数为：
该势函数说明，此MRF模型偏好你变量xA与xC拥有相同的取值，xB与xC拥有不同的取值；即，在该模型中xA与xC正相关，xB与xC负相关。结合HammersleyClifford 定理可知，令xA与xC相同，xB与xC不同的变量值指派将取得较高的联合概率。
为了满足非负性，定义势函数通常采用指数函数，即：
其中αuv和βv是参数。系数αuv对应的项考虑了每一对结点的关系；系数βv对应的项仅考虑单结点。
1.8，code：
图像分割前后的效果比较：
# 基于马尔科夫随机场的图像分割：
import cv2
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
%matplotlib inline
import warnings
warnings.filterwarnings("ignore") # 屏蔽警告
img = cv2.imread('../DataSets/ZiLong.jpg')
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) #将图片二值化，彩色图片该方法无法做分割
img = gray
img_double = np.array(img, dtype = np.float64)
cluster_num = 2
max_iter =200
label=np.random.randint(1, cluster_num + 1, size = img_double.shape)
iter=0
f_u = np.array([0,1,0,0,0,0,0,0,0]).reshape(3, 3)
f_d = np.array([0,0,0,0,0,0,0,1,0]).reshape(3, 3)
f_l = np.array([0,0,0,1,0,0,0,0,0]).reshape(3, 3)
f_r = np.array([0,0,0,0,0,1,0,0,0]).reshape(3, 3)
f_ul = np.array([1,0,0,0,0,0,0,0,0]).reshape(3, 3)
f_ur = np.array([0,0,1,0,0,0,0,0,0]).reshape(3, 3)
f_dl = np.array([0,0,0,0,0,0,1,0,0]).reshape(3, 3)
f_dr = np.array([0,0,0,0,0,0,0,0,1]).reshape(3, 3)
while iter < max_iter:
iter = iter + 1
#print(iter)
label_u = cv2.filter2D(np.array(label, dtype = np.uint8), -1 , f_u)
label_d = cv2.filter2D(np.array(label, dtype = np.uint8), -1 , f_d)
label_l = cv2.filter2D(np.array(label, dtype = np.uint8), -1 , f_l)
label_r = cv2.filter2D(np.array(label, dtype = np.uint8), -1 , f_r)
label_ul = cv2.filter2D(np.array(label, dtype = np.uint8), -1 , f_ul)
label_ur = cv2.filter2D(np.array(label, dtype = np.uint8), -1 , f_ur)
label_dl = cv2.filter2D(np.array(label, dtype = np.uint8), -1 , f_dl)
label_dr = cv2.filter2D(np.array(label, dtype = np.uint8), -1 , f_dr)
m, n = label.shape
p_c = np.zeros((cluster_num, m , n))
for i in range(cluster_num):
label_i = (i+1) * np.ones((m, n))
u_T = 1 * np.logical_not(label_i - label_u)
d_T = 1 * np.logical_not(label_i - label_d)
l_T = 1 * np.logical_not(label_i - label_l)
r_T = 1 * np.logical_not(label_i - label_r)
ul_T = 1 * np.logical_not(label_i - label_ul)
ur_T = 1 * np.logical_not(label_i - label_ur)
dl_T = 1 * np.logical_not(label_i - label_dl)
dr_T = 1 * np.logical_not(label_i - label_dr)
temp = u_T + d_T + l_T + r_T + ul_T + ur_T + dl_T + dr_T
p_c[i, :] = (1.0/8) * temp
p_c[p_c == 0] = 0.001
mu = np.zeros((1, cluster_num))
sigma = np.zeros((1, cluster_num))
for i in range(cluster_num):
index = np.where(label == (i+1))
data_c = img[index]
mu[0, i] = np.mean(data_c)
sigma[0, i] = np.var(data_c)
p_sc = np.zeros((cluster_num, m , n))
one_a = np.ones((m, n))
for j in range(cluster_num):
MU = mu[0, j]* one_a
p_sc[j, :] = (1.0/np.sqrt(2 * np.pi * sigma[0, j])) * np.exp(-1. * (( img - MU)**2) /(2 * sigma[0, j]))
X_out = np.log(p_c) + np.log(p_sc)
label_c = X_out.reshape(2, m * n)
label_c_t = label_c.T
label_m = np.argmax(label_c_t, axis = 1)
label_m = label_m + np.ones(label_m.shape) # 由于上一步返回的是index下标，与label其实就差1，因此加上一个ones矩阵即可
label= label_m.reshape(m, n)
label = label - np.ones(label.shape) # 为了出现0
lable_w = 255 * label #此处做法只能显示两类，一类用0表示另一类用255表示
cv2.imwrite('../OutPut/label_ZiLong.jpg', lable_w)
im_before = Image.open('../DataSets/ZiLong.jpg')
im_after = Image.open('../OutPut/label_ZiLong.jpg')
figsize = (16, 8)
f, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=figsize)
ax1.imshow(im_before)
ax2.imshow(im_after)
f.show()本文分享自 MiningAlgorithms 微信公众号，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体分享计划 ，欢迎热爱写作的你一起参与！