基于马尔可夫随机场(MRF)的遥感图像分割方法介绍
基于马尔可夫随机场(MRF)的遥感图像分割方法介绍
Astrophel背景与概述
在遥感影像的应用领域，图像分割是必不可少的步骤之一。传统的图像分割算法基于人工阈值，受主观影响大，造成效率低下；自动阈值法如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
62function 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 |
