基于马尔可夫随机场(MRF)的遥感图像分割方法介绍

背景与概述

在遥感影像的应用领域,图像分割是必不可少的步骤之一。传统的图像分割算法基于人工阈值,受主观影响大,造成效率低下;自动阈值法如OTSU,虽然解决了人工阈值确定法的缺点,但仍要求图像直方图近似双峰分布,对类别分布不均,感兴趣类别占比小的遥感影像分割效果不佳。
马尔可夫随机场(MRF)是一种基于贝叶斯理论的图像分割方法,该方法充分考虑了图像的邻域信息,且无需人工干预,具有效率高,精度好的优点。

马尔可夫随机场

原理

要理解马尔可夫随机场模型进行图像分割的原理,首先要认识到“马尔可夫性(Markov property)”。马尔可夫性是指“当一个随机过程在给定现在状态及所有过去状态情况下,其未来状态的条件概率分布仅依赖于当前状态”。用常见的天气作为例子:若某地明天的天气仅与今天的天气有关,而与昨天,前天都不相关,则认为该地天气具有“马尔可夫性”。
具体到图像分割,可以假设图像的某一像素所属类别只与其周围邻域的其他像素有关,而与邻域外的像素无关,即图像具有“马尔可夫性”。事实上包括遥感影像在内的大部分图像确实具有类似的特征,因此适合使用MRF模型进行图像分割。
邻域
假设待分割图像为S,分割结果为W,图像分割问题可以被转换为一个概率问题:即已知S中某一点的值,求解其在W中属于每一类的概率,取最大的概率所属类别作为该点的类别。例如,将一幅图像分割为两类,该图像中某一像素属于类别一的$P(W|S)=0.2$,属于类别二的$P(W|S)=0.4$,则该像素为类别二。那么,如何求解$P(W|S)$的最大值呢?
根据贝叶斯公式:

由于观测场(图像)已知,因此$P(S)$为定值,求解$P(W|S)$的最大值就是求解$P(S|W)P(W)$的最大值。
$P(S|W)$是指已知某像素属于W,其值为S的条件概率。假设类别概率分布情况已知,将S的值分别代入概率分布密度函数即可求出S在每一类中$P(S|W)$的值。
$P(W)$是该点的先验概率,该值理论上可用Hammcrslcy-Clifford定理求解,实际应用中,可以用迭代的方法逐渐逼近,即给予该点一个随机的先验概率,通过不断迭代,逐步得到接近正确的$P(W)$值。最后,将二者相乘,计算该点在每个个类别中的$P(S|W)P(W)$,取最大值所属于的类别作为该点的类别。

计算过程

下面将以一个像素及其邻域作为案例,计算一步MRF图像分割所进行的计算:
假设图像中某一像素及其八邻域如下矩阵所示:

初始标签场采用随机场,该点及其邻域的标签如下矩阵所示:

首先计算$P(W)$

然后计算$P(S|W)$,假设图像各类别服从高斯分布,且参数分别为$\mu_1 = 0.2, \sigma_1 = 0.1; \mu_2 = 0.8, \sigma_2 = 0.2$根据高斯分布的概率密度函数:

分别代入$x=0.9$,得:

此时,可求出:

因此,该点所属类别为2,更新后的标签阵如下矩阵所示:

依次更新其他位置的标签,然后迭代,逐步逼近最终结果。
这种情况较为理想,事实上,采用随机初始标签场的高斯分布参数差异不大,需要经过多次迭代才能出现明显的差异。因此,可以采取一些措施优化,减少迭代次数。

优化与改进

使用随机标签场(即随机$P(W)$)进行迭代求解,虽然较为方便,但迭代次数多,结果精度受限。因此,可以使用预分割的结果作为初始标签场进行迭代。相比于随机标签场,使用初始标签场的MRF具有可以减少迭代次数、控制输出类别的值等优势。常见的初始标签场可以为K-Means结果、OTSU结果等。

实现过程

使用Matlab实现上述计算过程。与分步计算案例一致:假设图像每一类的分布情况服从高斯分布,可以通过$\sigma$和$\mu$描述其分布,带入图像每一点的像素值就可以求出$P(S|W)$;$P(W)$的计算方法则采用八邻域法,通过计算每一个像素点周围的八个像素点属于每种类别的数量,求和除8,得到该点属于每个类别的概率。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
function bw = mrf2D(img,classnum,label,maxnum)
%mrf2D 使用马尔可夫随机场对图像进行分割
% img-被分割图像
% classnum-分割类别数
% label-初始标记场,可使用预分割结果,也可初始化一个随机场
% maxnum-最大迭代次数,避免过拟合
label = label + 1; % 输入label是从0开始
diff = 9999999;
num = 0;
while diff > 50000
% 计算p(w)
label_u = imfilter(label,[0,1,0;0,0,0;0,0,0],'replicate'); % 取8邻域值
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');
pw = zeros(classnum,size(label,1)*size(label,2)); % 初始化P(w)
for i = 1:classnum
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);
pw(i,:) = temp(:)/8; % 计算概率
end
% 计算P(s|w)
mu = zeros(1,classnum);
sigma = zeros(1,classnum);
%求出每一类的的高斯参数
for i = 1:classnum
index = label == i; %找到每一类的点
data_c = img(index);
data_c = data_c(~isnan(data_c)); % 去除数组中的NaN值
mu(i) = mean(data_c); %均值
sigma(i) = var(data_c);% 方差
end
p_sw = zeros(classnum,size(label,1)*size(label,2)); % 初始化P(s|w)
for i = 1:size(img,1)*size(img,2)
for j = 1:classnum
p_sw(j,i) = 1/sqrt(2*pi*sigma(j)) * exp(-(img(i)-mu(j))^2/2/sigma(j)); % 假设样本服从正态分布,带入概率密度函数求该点在每一类样本中的概率
end
end
p = pw .* p_sw; % 相乘得MRF模型贝叶斯公式的分子部分
% 求相乘最大值
label_1 = zeros(size(label));
for i = 1:size(img,1)*size(img,2)
if p(1,i) >= p(2,i)
label_1(i) = 1;
else
label_1(i) = 2;
end
end
% 迭代终止条件若干
diff = sum(sum(abs(label-label_1)));
label = label_1;
num = num + 1;
if num >= maxnum
break;
end
end
bw = label;
end

总结

基于MRF模型的图像分割方法有效利用了邻域信息,不要求图像直方图呈现明显的双峰分布,且可以自由切换其他概率分布模型,具有适应性强的优点。
经实验验证,MRF模型可以有效提升图像分割的正确率:

方法 OA Kappa
OTSU 89.29% 0.7836
MRF 90.12% 0.8021
OTSU-MRF 90.50% 0.8087

参考

  1. CSDN-马尔科夫随机场之图像分割【二】
  2. CSDN-从贝叶斯理论到马尔可夫随机场(MRF)—以图像分割为例