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 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!