MRF.zip_MATLAB 马尔可夫_mrf_图像分割 mrf_马尔可夫_马尔可夫 分割

function [imPara]=HmtTrain(I,L) %将图像I进行采样后，进行L层小波分解，计算参数imPara %imPara为结构体数组，每一元素对应每一尺度上的参数集合 %imPara（n）包括H、V、D。而每一个H、D、V又分别是一个结构体 %包括es,ps和si三个元素 I=double(I); %采样图像I，生成10个采用图像 sample=GetSample(I,10,64); %计算采样图像的小波变换 w=repmat(struct('wn',[]),[10,1]); for n=1:10 In=sample(n).data; w(n).wn=MyHaar(In,L); end %初始化 imPara=Initialization(w,L); %HMT训练 ps=struct('P1',[],'P2',[]); P=repmat(ps,[10,L]); %水平方向频域分量训练 maxiter=300;epsilon=1e-5; iter=0;err=Inf; error=zeros(maxiter,1); while (iter<maxiter)&&(err>epsilon) P=Expectation_H(w,imPara,P); [imParaN]=Maximization_H(w,P,imPara); err=error_H(imPara,imParaN); imPara=imParaN; iter=iter+1; error(iter)=err; end %垂直方向频域分量训练 iter=0;err=Inf; while (iter<maxiter)&&(err>epsilon) P=Expectation_V(w,imPara,P); [imParaN]=Maximization_V(w,P,imPara); err=error_V(imPara,imParaN); imPara=imParaN; iter=iter+1; error(iter)=err; end %对角线方向频域分量训练 iter=0;err=Inf; while (iter<maxiter)&&(err>epsilon) P=Expectation_D(w,imPara,P); [imParaN]=Maximization_D(w,P,imPara); err=error_D(imPara,imParaN); imPara=imParaN; iter=iter+1; error(iter)=err; end end % GetSample(I,sampleNumber,winsize)函数用于在打开的图像I中随机选取sampleNumber个大小为winsize的样本图像 function sample=GetSample(I,sampleNumber,winsize) %获取采样数据 %从图像I中随机获取sampleNumber个大小为winsize*winsize的样本数据，存入sample中 sizex=size(I,2);sizey=size(I,1); posx=floor(rand(1,sampleNumber)*(sizex-winsize-2))+1; posy=floor(rand(1,sampleNumber)*(sizey-winsize-2))+1; sample=struct('data',[]); sample=repmat(sample,[sampleNumber,1]); for j=1:sampleNumber sample(j).data=I(posy(1,j):posy(1,j)+winsize-1,posx(1,j):posx(1,j)+winsize-1,:); end end %函数MyHaar(I,L)用于计算图像I的L层Haar小波变换，并将变换结果用结构体数组进行存储 function w=MyHaar(I,L) w=repmat(struct('H',[],'V',[],'D',[]),[L,1]); for n=1:L [A,H,V,D]=dwt2(I,'Haar'); w(n).H=H; w(n).V=V; w(n).D=D; I=A; end end function [imPara]=Initialization(w,L) % 初始化：根据w中的小波系数取值，计算es,ps,si的初始取值 parab=struct('es',[],'ps',[],'si',[]); paras.H=parab;paras.V=parab;paras.D=parab; imPara=repmat(paras,[L,1]); for s=1:L es=ones(2)/2; ps=ones(2,1)/2; imPara(s).H.es=es;imPara(s).V.es=es;imPara(s).D.es=es; imPara(s).H.ps=ps;imPara(s).V.ps=ps;imPara(s).D.ps=ps; sih=0;siv=0;sid=0; for n=1:10 wn=w(n).wn; H=wn(s).H;V=wn(s).V;D=wn(s).D; sih=sih+mean(H(:).^2); siv=siv+mean(V(:).^2); sid=sid+mean(D(:).^2); end sih=sih/10;siv=siv/10;sid=sid/10; sih=sih*(sih>1e-6)+1e-6*(sih<=1e-6); siv=siv*(siv>1e-6)+1e-6*(siv<=1e-6); sid=sid*(sid>1e-6)+1e-6*(sid<=1e-6); imPara(s).H.si=sih*linspace(0.5,2.0,2)'; imPara(s).V.si=siv*linspace(0.5,2.0,2)'; imPara(s).D.si=sid*linspace(0.5,2.0,2)'; end end function [P]=Expectation_H(w,imPara,P) L=size(w(1).wn,1); for n=1:10 wn=w(n).wn; P(n,:)=E_H(wn,imPara,P(n,:)); end end function [P]=E_H(w,imPara,P) %计算小波系数w的H部分的期望 L=size(w,1); %定义中间变量 belt=repmat(struct('beltn',[]),[L,1]); beltip=repmat(struct('beltipn',[]),[L,1]); beltpi=repmat(struct('beltpin',[]),[L,1]); alpha=repmat(struct('alphan',[]),[L,1]); %最高分辨率尺度的初始化 n=1; wn=w(n).H; es=imPara(n).H.es;ps=imPara(n).H.ps;si=imPara(n).H.si; belt(n).beltn=gauss(wn,si); % up过程 for n=1:L-1 wn=w(n).H; [Widthn,Heightn]=size(wn); es=imPara(n).H.es;ps=imPara(n).H.ps;si=imPara(n).H.si; wp=w(n+1).H; [Widthp,Heightp]=size(wp); % 计算beltip beltn=belt(n).beltn; beltn=reshape(beltn,[2,Widthn*Heightn]); beltipn=es'*beltn; beltipn=reshape(beltipn,[2,Widthn,Heightn]); beltip(n).beltipn=beltipn; %计算belt(n+1) sp=imPara(n+1).H.si; gausp=gauss(wp,sp); tm=ones(2,Widthp,Heightp); tm=tm.*beltipn(:,1:2:Widthn,1:2:Heightn); tm=tm.*beltipn(:,2:2:Widthn,1:2:Heightn); tm=tm.*beltipn(:,1:2:Widthn,2:2:Heightn); tm=tm.*beltipn(:,2:2:Widthn,2:2:Heightn); beltp=gausp.*tm; belt(n+1).beltn=beltp; %计算beltpi tm=zeros(2,Widthn,Heightn); tm(:,1:2:Widthn,1:2:Heightn)=beltp; tm(:,2:2:Widthn,1:2:Heightn)=beltp; tm(:,1:2:Widthn,2:2:Heightn)=beltp; tm(:,2:2:Widthn,2:2:Heightn)=beltp; beltpi(n).beltpin=tm./beltipn; end %最高尺度初始化 n=L; wn=w(n).H; [Widthn,Heightn]=size(wn); ps=imPara(n).H.ps; alpha(n).alphan=repmat(ps,[1,Widthn,Heightn]); % down过程 for n=L-1:-1:1 wn=w(n).H; [Widthn,Heightn]=size(wn); es=imPara(n).H.es; %计算alphan beltpin=beltpi(n).beltpin; beltpin=reshape(beltpin,[2,Widthn*Heightn]); alphap=alpha(n+1).alphan; tm=zeros(2,Widthn,Heightn); tm(:,1:2:Widthn,1:2:Heightn)=alphap; tm(:,2:2:Widthn,1:2:Heightn)=alphap; tm(:,1:2:Widthn,2:2:Heightn)=alphap; tm(:,2:2:Widthn,2:2:Heightn)=alphap; alphap=tm; alphap=reshape(alphap,[2,Widthn*Heightn]); tm=es*(alphap.*beltpin); alpha(n).alphan=reshape(tm,[2,Widthn,Heightn]); end %计算P for n=1:L-1 wn=w(n).H; [Widthn,Heightn]=size(wn); es=imPara(n).H.es; %计算P1 alphan=alpha(n).alphan; beltn=belt(n).beltn; tm=repmat(sum(alphan.*beltn,1),[2,1,1]); P(n).P1=alphan.*beltn./tm; %计算P2 alphap=alpha(n+1).alphan; tm1=zeros(2,Widthn,Heightn); tm1(:,1:2:Widthn,1:2:Heightn)=alphap; tm1(:,2:2:Widthn,1:2:Heightn)=alphap; tm1(:,1:2:Widthn,2:2:Heightn)=alphap; tm1(:,2:2:Widthn,2:2:Heightn)=alphap; alphap=tm1; alphap=reshape(alphap,[2,Widthn*Heightn]); beltpin=beltpi(n).beltpin; beltpin=reshape(beltpin,[2,Widthn*Heightn]); for s=1:2 for p=1:2 belts=beltn(s,:); esp=es(s,p); alphaps=alphap(p,:); beltpips=beltpin(p,:); fenmu=tm(s,:); tmp=belts*esp.*alphaps.*beltpips./fenmu; tmp=reshape(tmp,[Widthn,Heightn]); P(n).P2(s,p,:,:)=tmp; end end end n=L; alphan=alpha(n).alphan; beltn=belt(n).beltn; tm=repmat(sum(alphan.*beltn,1),[2,1,1]); P(n).P1=alphan.*beltn./tm; end function [y]=gauss(w,s) y=[]; for m=1:2 y=cat(3,y,1./sqrt(2.0*pi*s(m)).*exp(-w.^2./(2.0*s(m)))+1e-15); end y=shiftdim(y,2); end function [P]=Expectation_V(w,imPara,P) L=size(w(1).wn,1); for n=1:10 wn=w(n).wn; P(n,:)=E_V(wn,imPara,P(n,:)); end end function [P]=E_V(w,imPara,P) %计算小波系数w的V部分的期望 L=size(w,1); %定义中间变量 belt=repmat(struct('beltn',[]),[L,1]); beltip=repmat(struct('beltipn',[]),[L,1]); beltpi=repmat(struct('beltpin',[]),[L,1]); alpha=repmat(struct('alphan',[]),[L,1]); %最高分辨率尺度的初始化 n=1; wn=w(n).V; es=imPara(n).V.es;ps=imPara(n).V.ps;si=imPara(n).V.si; belt(n).beltn=gauss(wn,si); % up过程 for n=1:L-1 wn=w(n).V; [Widthn,Heightn]=size(wn); es=imPara(n).V.es;ps=imPara(n).V.ps;si=imPara(n).V.si; wp=w(n+1).V; [Widthp,Heightp]=size(wp); % 计算beltip beltn=belt(n).beltn; beltn=reshape(beltn,[2,Widthn*Heightn]); beltipn=es'*beltn; beltipn=reshape(beltipn,[2,Widthn,Heightn]); beltip(n).beltipn=beltipn; %计算belt(n+1) sp=imPara(n+1).V.si; gausp=gauss(wp,sp); tm=ones(2,Widthp,Heightp); tm=tm.*beltipn(:,1:2:Widthn,1:2:Heightn); tm=tm.*beltipn(:,2:2:Widthn,1:2:Heightn); tm=tm.*beltipn(:,1:2:Widthn,2:2:Heightn); tm=tm.*beltipn(:,2:2:Widthn,2:2:Heightn); beltp=gausp.*tm; belt(n+1).beltn=beltp; %计算beltpi tm=zeros(2,Widthn,Heightn); tm(:,1:2:Widthn,1:2:Heightn)=beltp; tm(:,2:2:Widthn,1:2:Heightn)=beltp; tm(:,1:2:Widthn