从贝叶斯理论到马尔可夫随机场--以图像分割为例

马尔可夫随机场（CRF）图像分割过程Matlab代码实现Python实现代码参考文献

本文主要介绍马尔可夫随机场及其在图像分割中的应用。基于马尔可夫随机场（MRF）的图像分割是一种基于统计学习方法的图像分割算法，其模型参数较少，空间约束性较强，使用较广泛。

马尔可夫随机场（CRF）

首先了解以下马尔可夫模型，纯粹的马尔可夫模型是指一个事物的当前状态只与之前的1个状态有关，而与其他状态没有关系。比如今天的天气好坏只与昨天的天气有关，而与前天以及大前天都没有关系，符合这样的一种特性的事物认为其具有马尔可夫性。

那么引申到图像领域，就是认为图像中某一点的特征（一般都是像素点的灰度值或颜色值）只与其附近的一小块邻域有关，而与其他区域无关。

图像分割过程

马尔可夫随机场在图像领域的重要应用就是图像分割。关于图像分割问题，从聚类角度讲就是一个图像聚类问题，把具有相同性质的像素点聚集为一类。简而言之，图像分类就是像素级别的分类问题，比如把一个图像分割为四个部分，就是把所有的像素点分为四类，即给每个像素点找一个标签类。

假设待分割的图像是S，大小为m*n，那么把每个像素点放到集合S中，图像就是S={S1,S2,⋯,Sm∗n}S=\{S_1,S_2,\cdots,S_{m*n}\}S={S1​,S2​,⋯,Sm∗n​}，把待分割的图像称为观测图像。假设图像已经分割完成，每个像素都有具体的类别，把最终的分割结果称为W，显然W的shape与S的shape相同，每个类别的取值范围是1-4.

现在的问题就是如何求W，我们所知道就是观测到的图像S，也就是说知道S的情况下求W，转化为概率就是求取P(W∣S)P(W|S)P(W∣S)，图像分割问题就变成了求取这个概率的最大值，也就是说根据S计算出图像的分割标签。

根据贝叶斯公式P(Y∣X)=P(X∣Y)P(Y)P(X)P(Y|X)={{P(X|Y)P(Y)}\over{P(X)}}P(Y∣X)=P(X)P(X∣Y)P(Y)​.

W为要求的分类标签，P(W)为标记场w的先验概率，P(S|W)是观察值S的条件概率分布，也称为似然函数。在已知像素点标记w的情况下，真实的像素点的概率，所以是一个似然函数，所以P(S)认为是个定值。那么现在要求P(W|S)的最大值，也就是要求P(S|W)P(W)的最大值。

先来说说P(W)这个标记场w的先验概率。那么我们落实到图像的每一个像素点上，也就是说这个像素点是标签L的概率是多少，有人可能会说，分割之前图像的分类标签都不知道，还怎么算是某一类标签L的概率，这个问题是有点较劲不好理解。

我觉得，这就是一个鸡蛋问题，是先有鸡还是先有蛋的问题，我们要求这只鸡，但是又没有蛋怎么办？一个简单的办法是我随便弄一个蛋来，管他是鸡蛋鸭蛋，有了蛋，鸡（或者鸭）就可以有了，虽然这个鸡不太像鸡，比如说是只野鸡，那么有了野鸡，野鸡再稍微进化一下，然后再下个蛋，蛋又长出野鸡1号，野鸡1号再稍微进化一下，然后再下个蛋，变成野鸡2号，如此反复反复，知道野鸡n号，然后蓦然回首发现，野鸡n号和我们所要的鸡竟是那么的相似，那么把它就认为是我们要的鸡了。有没有糊涂？

其实优化聚类里面很多问题都是鸡蛋问题，比如曾经介绍的FCM算法，有个先给u还是先给c的鸡蛋问题，u的计算需要c，c的计算有需要u，那怎么办，先假设一个吧。再比如EM算法的迭代过程也可以看成鸡蛋问题，有了E步算M步，在回来算E步，在M步。。。。最终得到最优解。

我们的问题要干嘛？要求一个像素点是标签L的概率，开始标签没有，ok假设一个吧，一个不行，那么在开始我们把整幅图像初始化一个分割标签，尽管这个标签不对（野鸡），也可以假设一个稍微对一点的标签（比如用其他分类算法（kmeans）得到一个预分割标签）（这个时候认为是野鸡100号）（好处在于这个时候算法不用迭代那么多次就可以停止了）。好了有了初始的标签，那么再来求一个像素点是标签L的概率。从这个时候就需要马尔科夫协助了。

我们想，要求一个像素点是标签L的概率，此时这个像素点附近的领域都也已经划分标签了，虽然这个像素点也有一个标签，但是这个标签是上一代的标签，那么它到下一代是需要进化的，马尔科夫性质告诉我们这个像素点的分类情况只与附近一些领域分类情况有关，而与另一些领域没有关系，也就是说我们可以根据这个像素点的附近领域的分类情况决定这个像素点在下一代是属于哪一类的。

既然我们认为每个像素点的分类符合马尔科夫随机模型，而另一些学者证明了这种马尔科夫随机场可以与一个Gibbs随机场等价（这就是Hammcrslcy-Clifford定理，人家证明的，那就这样叫吧）。而Gibbs随机场也有一个概率密度函数，这样我们就用求图像的Gibbs随机场的概率P代替了P(W)。那么Gibbs随机场的概率P怎么表示呢？如下：

P(w)=1zexp(−1TU2(w))P(w)={1\over z}exp(-{1\over T}U_2(w))P(w)=z1​exp(−T1​U2​(w))

关于P(W)就说这么多，下面说说P(S|W)。P(S|W)已知分类标签那么它的像素值（灰度）是s的概率。现在就假设w=1，某个像素点灰度为s，那么表示的意思就是在第一类里面像素灰度为s的概率。因为分类标签在前面说到，每次迭代的时候是有一个分类标签的（尽管不是最后的标签），那么我们可以把属于第一类的所有点都挑出来，考虑每个点都是独立的，并且认为每一类里面的所有点服从高斯分布（正态分布），那么在每一类里面我们是不是可以根据这一类里面的这些点建立一个属于这一类的高斯密度函数了，那么再来一个像素点值，把它带到这类里面密度函数中去就可以得到这个概率了吧。同理对于2,3,4类，每一类都可以建立一个高斯密度函数，这样就有四个高斯密度函数了，那么每一个点属于这四类的概率就可以分别带到这四个高斯密度函数中计算了。高斯密度函数一般形式为：

P(x∣wi)=12πσexp(−(x−u)22σ2)P(x|w_i)={1\over{\sqrt{2\pi}\sigma}}exp(-{{(x-u)^2}\over{2\sigma^2}})P(x∣wi​)=2π​σ1​exp(−2σ2(x−u)2​)

举个例子假设现在我们求得的四类高斯函数，其图像如下所示：

某一点的灰度为s=110，那么它对应的分别P(s|W1)、P(s|W2)、P(s|W3)、P(s|W4)，就分别如上图所示呢，很显然的可以看到110在第二类下的概率最大，其他的几个概率也可以出来的。

现在这一部分对于每个点也有了四个概率，结合上面每个点的四个概率，那么对每个点，属于每一类的概率对应相乘，得到最终每个点属于四个类的概率。这个时候就可以挑其中最大的那么就是该点所属于的类了。

这样一次迭代，所有的点所属于的类更新一遍，那这个新的类标签作为下一次迭代的类标签，反复下去，程序结束的条件可以是设置迭代次数，也可以是观察类中心不在变化为止。

好了，上述的概率相乘取最大是原始一点的算法。其实可以看到，这个计算包括两个部分，一般的讲法，认为这两部分每一部分都组成一个能量，换个说法就是最优化能量函数了，可以表示为如下：

W=argmin(U1(w,S),U2(w))W=arg\space min(U_1(w,S),U_2(w))W=argmin(U1​(w,S),U2​(w))

像上述的概率相乘在实际中取一下对数就变成相加了。更深层次的马尔科夫随机场可以去看看相关论文，虽然方法可能不太一样，但是原理上是一样的。

Matlab代码实现

下面以条件迭代算法（ICM）来实例阐述MRF在图像分割上的应用，编程平台为matlab，相关代码如下：

首先将彩色图像转换为灰度图像

img_rgb=imread('lena.jpg'); %读入需要转换的彩色图像imshow(img_rgb); %显示彩色图像img_gray=rgb2gray(img_rgb); %RGB图像转为灰度图像figureimshow(img_gray); %显示灰度图像imwrite(img_gray,'lena_gray.jpg'); %保存灰度图像

使用K-Means对图像进行分割，类别数为2

clcclearclose allimg = double(imread('lena_gray.jpg')); %more quicklycluster_num = 2; %设置分类数%-----------kmeans初始化预分割----------label = kmeans(img(:),cluster_num);label = reshape(label,size(img));imshow(label,[]);title('kmeans初始化预分割');

基于马尔可夫随机场的图像分割代码

clcclearclose allimg = double(imread('lena_gray.jpg')); %more quicklycluster_num = 4; %设置分类数maxiter = 60; %最大迭代次数%-------------随机初始化标签----------------label = randi([1,cluster_num],size(img));%-----------kmeans最为初始化预分割----------% label = kmeans(img(:),cluster_num);% label = reshape(label,size(img));iter = 0;while iter < maxiter%-------计算先验概率---------------%这里我采用的是像素点和3*3领域的标签相同%与否来作为计算概率%------收集上下左右斜等八个方向的标签--------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领域标签相对于每一类的相同个数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;%防止出现0%---------------计算似然函数----------------mu = zeros(1,4);sigma = zeros(1,4);%求出每一类的的高斯参数--均值方差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%------计算每个像素点属于每一类的似然概率--------%------为了加速运算，将循环改为矩阵一起操作--------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 %找到联合一起的最大概率最为标签，取对数防止值太小[~,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

当采用随机初始化分割标签，类别数为4类，迭代次数为60次时的分割结果如下。

当采用K-Means初始化分割标签，类别数为4类，迭代次数为60次时的分割结果如下。

当采用K-Means初始化分割标签，类别数为2类，迭代次数为60次时的分割结果如下。

上述程序除了注释外再说一点，那就是对于是否需要预分类，程序里面写了随机分类和kmeans预分类两者。当分别去实验可以发现，如果采用kmeans预分类后，迭代开始分割的结果就很好，程序会在这个基础上再变化一点。采用kmeans预分类的好处在于，分割的结果会更快达到稳定，且更准确。

而采用随机的预分类，分割的结果也算可以，这种方法得到的是一个局部最优解，当然kmeans得到的也是个局部最优解（不敢说最优解），但是前面的局部最优解比kmeans后的局部最优解又会差一点，尤其是在分割类别数大一点的时候，有kmeans预分类的效果会明显好于随机预分类。因为类别数越是大，相当于维度就越是大，那么它存在的局部最优解就越是多，那么从随机的预分类（最差的情况）往最优解方向走时，中途可能随便碰到一个局部最优解就走不动了。从某种意义上讲，这个分割是一个np难问题，很难找到分割的最优解。

Python实现代码

将彩色图片转为灰度图像

import cv2import 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()

参考文献

/on2way/article/details/47307927