Matlab程序代码
function [seg]=MRMRF(w,class_number,potential,maxIter)
%MRMRF图像分割算法
%w-待分解图像的多尺度序列
%class_number-分类数
%potential-potts模型势函数
%maxIter-最大迭代次数
%多分辨率表达的尺度数
L=size(w,1);
%分割结果的多尺度序列
seg=cell(L,1);
%使用ICM算法计算最高尺度的分割结果
seg{
L
}=ICM(w{
L
},class_number,potential,maxIter);
%计算其他尺度的分割结果
for n=(L-1):-1:1
%获取初始分割结果
segn=myZoomOut(seg{
n+1
});
%获取尺度n上的最终分割结果
wn=w{
n
};
segn=ICMn(wn,segn,class_number,potential,maxIter);
%保存尺度n上的分割结果
seg{
n
}=segn;
end
end
function [seg]=ICMn(image,seg,class_number,potential,maxIter) %尺度n上已知初始结果的ICM分割
[width,height,bands]=size(image);
%将图像和初始分割结果分别转换为向量
image=imstack2vectors(image);
seg=imstack2vectors(seg);
%ICM迭代
iter=0;
while(iter<maxIter)
%估计特征场参数
[mu,sigma]=GMM_parameter(image,seg,class_number);
%计算特征场能量
Ef=EnergyOfFeatureField(image,mu,sigma,class_number);
%计算标记场能量
El=EnergyOfLabelField(seg,potential,width,height,class_number);
% 计算新的分割结果
E=Ef+El;
[tm,seg]=min(E,[],2);
% 更新迭代条件
iter=iter+1;
end
seg=reshape(seg,[width height]);
end
function [x]=myZoomOut(x)
%放大x两倍
[col,raw,B]=size(x);
tm=zeros(2*col,2*raw,B);
tm(1:2:2*col,1:2:2*raw,:)=x;
tm(2:2:2*col,1:2:2*raw,:)=x;
tm(1:2:2*col,2:2:2*raw,:)=x;
tm(2:2:2*col,2:2:2*raw,:)=x;
x=tm;
end
%高斯模型参数估计
function [mu,sigma]=GMM_parameter(image,seg,class_number)
% image-已经向量化的图像数据，每一行表示一个像素
% seg-与image对应的图像标记，每一个对应一个像素标记
% class_number-图像的分类数
[n,d]=size(image);
mu=zeros(class_number,d);
sigma=zeros(d,d,class_number);
for i=1:class_number
Im_i=image(seg==i,:);
[sigma(:,:,i),mu(i,:)]=covmatrix(Im_i);
end
end
function [C,m]=covmatrix(X)
[k,n]=size(X);
X=double(X);
if k==1
C=0;
m=X;
else
m=sum(X,1)/k;
X=X-m(ones(k,1),:);
C=(X'*X)/(k-1);
m=m';
end
end
function [E]=EnergyOfLabelField(seg,potential,width,height,class_number) % seg-向量化的像素标记矩阵
% potential-Potts模型势函数取值
% width,height-图像大小
n=size(seg,1);
seg=reshape(seg,[width height]);
%计算领域
nei8_image=neighbor8syn(seg);
Nei8=imstack2vectors(nei8_image);
%计算标记场能量
E=zeros(n,class_number);
for i=1:class_number
E(:,i)=sum(Nei8~=i,2);
end
E=E*potential;
end
function nei=neighbor8syn(Y)
[s,t,P]=size(Y);
Yul=zeros(s,t,P);
Yul(2:s,2:t,:)=Y(1:s-1,1:t-1,:);
Yu=zeros(s,t,P);
Yu(2:s,:,:)=Y(1:s-1,:,:);
Yur=zeros(s,t,P);
Yur(2:s,1:t-1,:)=Y(1:s-1,2:t,:);
Yr=zeros(s,t,P);
Yr(:,1:t-1,:)=Y(:,2:t,:);
Ydr=zeros(s,t,P);
Ydr(1:s-1,1:t-1,:)=Y(2:s,2:t,:);
Yd=zeros(s,t,P);
Yd(1:s-1,:,:)=Y(2:s,:,:);
Ydl=zeros(s,t,P);
Ydl(1:s-1,2:t,:)=Y(2:s,1:t-1,:);
Yl=zeros(s,t,P);
Yl(:,2:t,:)=Y(:,1:t-1,:);
nei=cat(3,Yul,Yu,Yur,Yr,Ydr,Yd,Ydl,Yl);
end
function [E]=EnergyOfFeatureField(image,mu,sigma,class_number) % mu-均值矩阵
% sigma-方差矩阵
n=size(image,1);
E=zeros(n,class_number);
% 计算能量
for i=1:class_number
mu_i=mu(i,:);
sigma_i=sigma(:,:,i);
diff_i=image-repmat(mu_i,[n,1]);
E(:,i)=sum(diff_i*inv(sigma_i).*diff_i,2)+log(det(sigma_i)); end
end。