Matlab程序代码.docx

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