Matlab程序代码

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);

%计算其他尺度的分割结果

forn=(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

%估计特征场参数

[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);

fori=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);

ifk==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);

fori=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);

% 计算能量

fori=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