# 从贝叶斯理论到图像马尔科夫随机场

P(W|S)=P(S|W)P(W)P(S)

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

P(w)=z1e1TU2(w)

Vc(wc)=Vs,r(ws,wr)={ββws=wtwswt
β$\beta$为耦合系数，通常0.5-1。

211134212

P(x|wi)=12πσexp((xu)22σ2)

W=argmin(U1(w,S)+U2(w))

clc
clear
close all
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