首先解释几个基本概念:
先验概率: 基于历史经验得到的当前事件发生的概率
后验概率: 根据数据或证据得到的事件发生的概率 (由因到果)
似然概率: 已知事件发生的概率下数据或证据发生的概率 (由果到因)
————————————————————————————————————————————————
图像分割和目标识别都可以看成一个图像中各个像素点做label的过程。例如进行边缘检测,就是以一定阈值判断这些像素值属于标签集{edge,nonedge}中的哪一类。图形分割从概率的角度来说就是一种求各像素点属于每个label的概率,取其中最大的那个label的过程。
P(w|s)表示基于像素信息s的条件下属于一个标签的概率,也叫后验概率。因此,该过程就是最大化后验概率来寻找最优标签w。
P(s|w)就是似然概率,在已知像素点标记w的情况下,它是真实的像素点s的概率; P(w)是先验概率,W就是我们要求的分类标签, P(s)是一个常量。由于图像中某些区域的相邻像素点见会有相似的颜色,纹理之类的,那么马尔科夫随机场这种概率模型就能够得到很好的应用。基于马尔科夫的随机场(MRF)的图像分割是一种基于统计的图像分割算法。
首先了解一下马尔科夫模型,纯粹的马尔科夫模型就是指一件事物的当前状态只与它之前的1个或者n个状态有关,而与再之前的状态没有关系,比如今天天气好坏只与昨天天气有关,而与前天乃至大前天都没有关系。符合这样的一种特性的事物认为其具有马尔科夫性。那么引申到图像领域,就是认为图像中某一点的特征(一般都是像素点灰色、颜色值等)只与其附近的一小块领域有关,而与其他的领域无关。想想也是这样,大多数时候,图像中某一点像素和附近的像素是相关的,附近领域像素是黑的,那它八成也是黑的,很好理解。
关于图像分割问题,从聚类角度讲,就是一个图像的聚类问题,把具有相同性质的像素点设置为一类。在说详细点就是一个标签分类问题,比如把一副图像分割成4类,那么每一个像素点必定属于这四类中的某一类,假设四类为1,2,3,4类,那么分割就是给每个像素点找一个标签类。好了,假设我们的待分割图像是S,大小m*n,那么把每个像素点放到集合S中,图像就是:S=S1,S2,...,Sm∗n。把待分割的图像称为观测到的图像。
现在假设分割好了,每个像素点都分了类,那么把最终的分割结果称为W,那么显然W大小与S一样大,W=w1,w2,...,wm∗n,其中所有的w取值都在1-4之间,4是最大的分类数,每个w都是在1-4之间取值的。
这是我们假设知道了最终分割结果W,但是W在开始是不知道的?现在的问题就是如何求W,我们所知道的不过是观测到的图像S,也就是说在知道S情况下求W,转化为概率就是求取P(W|S),问题就变成求取这个概率的最大值,也就是说根据S我们算出这个图像的最有可能的分割标签是什么。即:
先来说说P(W)这个标记场w的先验概率。那么我们落实到图像的每一个像素点上,也就是说这个像素点是标签L的概率是多少,有人可能会说,分割之前图像的分类标签都不知道,还怎么算是某一类标签L的概率,这个问题是有点较劲不好理解。我觉得,这就是一个鸡蛋问题,是先有鸡还是先有蛋的问题,我们要求这只鸡,但是又没有蛋怎么办?一个简单的办法是我随便弄一个蛋来,管他是鸡蛋鸭蛋,有了蛋,鸡(或者鸭)就可以有了,虽然这个鸡不太像鸡,比如说是只野鸡,那么有了野鸡,野鸡再稍微进化一下,然后再下个蛋,蛋又长出野鸡1号,野鸡1号再稍微进化一下,然后再下个蛋,变成野鸡2号,如此反复反复,知道野鸡n号,然后蓦然回首发现,野鸡n号和我们所要的鸡竟是那么的相似,那么把它就认为是我们要的鸡了。又没有糊涂?其实优化聚类里面很多问题都是鸡蛋问题。比如EM算法的迭代过程也可以看成鸡蛋问题,有了E步算M步,在回来算E步,在M步。。。。最终得到最优解。
回到我们的问题,我们的问题要干嘛?要求一个像素点是标签L的概率,开始标签没有,ok假设一个吧,一个不行,那么在开始我们把整幅图像初始化一个分割标签,尽管这个标签不对(野鸡),也可以假设一个稍微对一点的标签(比如用其他分类算法(kmeans)得到一个预分割标签)(这个时候认为是野鸡100号)(好处在于这个时候算法不用迭代那么多次就可以停止了)。好了有了初始的标签,那么再来求一个像素点是标签L的概率。从这个时候就需要马尔科夫协助了。我们想,要求一个像素点是标签L的概率,此时这个像素点附近的领域都也已经划分标签了,虽然这个像素点也有一个标签,但是这个标签是上一代的标签,那么它到下一代是需要进化的,马尔科夫性质告诉我们这个像素点的分类情况只与附近一些领域分类情况有关,而与另一些领域没有关系,也就是说我们可以根据这个像素点的附近领域的分类情况决定这个像素点在下一代是属于哪一类的。
Ok既然我们认为每个像素点的分类符合马尔科夫随机模型,而另一些学者证明了这种马尔科夫随机场可以与一个Gibbs随机场等价(这就是Hammcrslcy-Clifford定理,人家证明的,那就这样叫吧)。而Gibbs随机场也有一个概率密度函数,这样我们就用求图像的Gibbs随机场的概率P代替了P(W)。那么Gibbs随机场的概率P怎么表示呢?如下:
现在要计算中心的像素点在下一次迭代中是属于第几类(这一代是第3类),采用一阶势能,这里需要说明一点,这个像素点无非是1-4之间的一类,那么我们需要分别计算下一代它是第1,2,3,4类时的势团。如上图的3*3像素框,以中心为例,先假设这个点是第1类,先比较左、右、上、左下,发现都是1-1一样,记一下“β”,比较左上、右上、右下,发现与“1”不一样,那么记一下“−β”,这样就可以将中心点如果是第1类的势团计算出来。那么在计算中心点如果是第2类,发现这个时候除了3个斜着的方向是一样的外,其他的都不一样。在计算假设是第三类,第四类。这样有了每类下的势团,就可以计算概率了,最终得到这个点属于第一类的概率,第二类的概率,第三类的概率,第四类的概率,这四个值在后面会一起用到来决定这个点到底属于哪一类。你可能会说,知道了这四个概率,看哪个最大不就出来了吗?注意这还只是第一项的先验概率,后面还有个似然函数的概率没有计算呢,后面你会看到这个似然函数的概率对于每个像素点也有四个概率,在这四个概率都计算出来之后,那么这个点属于第一类的概率用着两个都属于第一类的概率相乘,第二类也是如此,最后再比较这四个概率的最大值作为要标记的标签。
通观Gibbs,其实就是一个求像素点与周围像素点的不一致程度,或者说这个像素点的周围像素点中,看看各个标签出现的概率多大,就是这个点属于对应标签的概率。所以在实际编程上,也没看到几个人实实在在的用它的原版公式编程,反而简化的计算,看看各个标签出现的次数等等方式来计算的。
下面说说P(S|W)。P(S|W),已知分类标签那么它的像素值(灰度)是s的概率。现在就假设w=1,某个像素点灰度为s,那么表示的意思就是在第一类里面像素灰度为s的概率。因为分类标签在前面说到,每次迭代的时候是有一个分类标签的(尽管不是最后的标签),那么我们可以把属于第一类的所有点都挑出来,考虑每个点都是独立的,并且认为每一类里面的所有点服从高斯分布(正态分布),那么在每一类里面我们是不是可以根据这一类里面的这些点建立一个属于这一类的高斯密度函数了,那么再来一个像素点值,把它带到这类里面密度函数中去就可以得到这个概率了吧。同理对于2,3,4类,每一类都可以建立一个高斯密度函数,这样就有四个高斯密度函数了,那么每一个点属于这四类的概率就可以分别带到这四个高斯密度函数中计算了。高斯密度函数一般形式为:
假设现在我们求得的四类高斯函数,其图像如下所示:
某一点的灰度为s=110,那么它对应的分别P(s|W1)、P(s|W2)、P(s|W3)、P(s|W3),就分别如上图所示呢,很显然的可以看到110在第三类下的概率最大,其他的几个概率也可以出来的。
现在这一部分对于每个点也有了四个概率,结合上面每个点的四个概率,那么对每个点,属于每一类的概率对应相乘,得到最终每个点属于四个类的概率。这个时候就可以挑其中最大的那么就是该点所属于的类了。
这样一次迭代,所有的点所属于的类更新一遍,那这个新的类标签作为下一次迭代的类标签,反复下去,程序结束的条件可以是设置迭代次数,也可以是观察类中心不在变化为止。
好了,上述的概率相乘取最大是原始一点的算法。其实可以看到,这个计算包括两个部分,一般的讲法,认为这两部分每一部分都组成一个能量,换个说法就是最优化能量函数了,可以表示为如下:
下面以条件迭代算法(ICM)来实例阐述MRF在图像分割上的应用,编程平台为matlab,相关代码如下:
clc
clear
close all
img = double(imread('lena.jpg'));%more quickly
cluster_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_num
label_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;%计算概率
end
p_c(find(p_c == 0)) = 0.001;%防止出现0
%---------------计算似然函数----------------
mu = zeros(1,4);
sigma = zeros(1,4);
%求出每一类的的高斯参数--均值方差
for i = 1:cluster_num
index = find(label == i);%找到每一类的点
data_c = img(index);
mu(i) = mean(data_c);%均值
sigma(i) = var(data_c);%方差
end
p_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_num
MU = 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;
end
subplot(2,3,n);
imshow(label,[])
title(['iter = ',num2str(iter)]);
pause(0.1);
n = n+1;
iter = iter + 1;
end
现在假设把图像分两类cluster_num=2,采用随机初始化分类,贴几个之间结果: