基于MRF多时相变化检测算法

tic t1=imread('After.png'); t2=imread('Before.png'); t1=imresize(t1,[328,338]); t2=imresize(t2,[328,338]); [a,b]=size(t1); [c,d]=size(t2); if size(t1,3)>1 t1=rgb2gray(t1); t2=rgb2gray(t2); end [hang,lie]=size(t1); r=(double(t1)+1)./(double(t2)+1);%求两幅图片的比值（只能检测灰度值变小的像素） Y=mat2gray(r); u=double(im2uint8(mat2gray(r)))+1;%比值图像 figure;imshow(Y,[]);title('比值图像'); T_yuzhi=0.35;%zuixiaocuofen(Y)/255;%最小错分概率求阈值 bianhua=Y>T_yuzhi;%求初次变化图像 figure;imshow(bianhua,[]);title('正变化初始分割'); window0=(bianhua==0); window1=(bianhua==1); k10=sum(sum(log(u).*window0))/sum(sum(window0));%不变部分均值,sum(sum(window0) 元素值的总和 k11=sum(sum(log(u).*window1))/sum(sum(window1));%变化部分均值 k20=sum(sum(((log(u)-k10).^2).*window0))/(sum(sum(window0))-1);%求不变部分的方差 k21=sum(sum(((log(u)-k11).^2).*window1))/(sum(sum(window1))-1);%求变化部分的方差 XX1=wextend('2D','sym',bianhua,[1,1]);%wextend 扩展矢量或矩阵，扩展1行1列 m0=zeros(hang,lie); m1=zeros(hang,lie); %%%%%%%%%%%%%%%%%%%%%%%%求k像素邻域相同标记的个数%%%%%%%%%%%%%% for i=2:hang+1 for j=2:lie+1 window3=XX1(i-1:i+1,j-1:j+1); if XX1(i,j)==0 m0(i-1,j-1)=8-sum(sum(window3));%8领域0的个数 m1(i-1,j-1)=sum(sum(window3));%8领域1的个数 else m0(i-1,j-1)=9-sum(sum(window3));%8领域0的个数 m1(i-1,j-1)=sum(sum(window3))-1;%8领域1的个数 end end end U0=((log(u)-k10).^2-k20*log(k20))/(2*k20)-0.3*m0;%求不变假设下的能量函数 U1=((log(u)-k11).^2-k21*log(k21))/(2*k21)-0.3*m1;%求变化假设下的能量函数 bianhua1=(U0>U1);%第一次得到的变化图像 chazhi=1; while chazhi>3e-3 p0=exp(-U0)./(exp(-U0)+exp(-U1));%不变部分的条件概率 p1=exp(-U1)./(exp(-U0)+exp(-U1));%变化部分的条件概率 w0=(p0>=p1).*p0; w1=(p1>p0).*p1; k10=sum(sum(w0.*log(u)))/(sum(sum(w0)));%求不变部分的均值 k11=sum(sum(w1.*log(u)))/(sum(sum(w1)));%求变化部分的均值 k20=sum(sum(((log(u)-k10).^2).*w0))/(sum(sum(w0)));%求不变部分的方差 k21=sum(sum(((log(u)-k11).^2).*w1))/(sum(sum(w1)));%求变化部分的方差 disp('求步长'); buchang=newton(w0,w1,m0,m1);%牛顿-拉夫逊算法求步长因子 XX1=wextend('2D','sym',bianhua1,[1,1]); for i=2:hang+1 for j=2:lie+1 window3=XX1(i-1:i+1,j-1:j+1); if XX1(i,j)==0 m0(i-1,j-1)=8-sum(sum(window3)); m1(i-1,j-1)=sum(sum(window3)); else m0(i-1,j-1)=9-sum(sum(window3)); m1(i-1,j-1)=sum(sum(window3))-1; end end end U0=((log(u)-k10).^2-k20*log(k20))/(2*k20)-buchang*m0; U1=((log(u)-k11).^2-k21*log(k21))/(2*k21)-buchang*m1; bianhua2=(U0>U1); chazhi=sum(sum(abs(bianhua1-bianhua2)))/(hang*lie); fprintf('连续两次的结果差值=%i\n',chazhi) bianhua1=bianhua2; % figure; % imshow(bianhua1); end figure;imshow(bianhua1,[]);title('正变化分割结果'); r=(double(t2)+1)./(double(t1)+1);%求两幅图片的比值（只能检测灰度值变大的像素） Y=mat2gray(r); u=double(im2uint8(mat2gray(r)))+1;%比值图像 % figure;imshow(Y,[]); T_yuzhi=zuixiaocuofen(Y)/255;%最小错分概率求阈值 bianhua=Y>T_yuzhi;%求初次变化图像 figure;imshow(bianhua,[]);title('反变化初始分割'); window0=(bianhua==0); window1=(bianhua==1); k10=sum(sum(log(u).*window0))/sum(sum(window0));%不变部分均值 k11=sum(sum(log(u).*window1))/sum(sum(window1));%变化部分均值 k20=sum(sum(((log(u)-k10).^2).*window0))/(sum(sum(window0))-1);%求不变部分的方差 k21=sum(sum(((log(u)-k11).^2).*window1))/(sum(sum(window1))-1);%求变化部分的方差 XX1=wextend('2D','sym',bianhua,[1,1]); m0=zeros(hang,lie); m1=zeros(hang,lie); %%%%%%%%%%%%%%%%%%%%%%%%求k像素邻域相同标记的个数%%%%%%%%%%%%%% for i=2:hang+1 for j=2:lie+1 window3=XX1(i-1:i+1,j-1:j+1); if XX1(i,j)==0 m0(i-1,j-1)=8-sum(sum(window3)); m1(i-1,j-1)=sum(sum(window3)); else m0(i-1,j-1)=9-sum(sum(window3)); m1(i-1,j-1)=sum(sum(window3))-1; end end end U0=((log(u)-k10).^2-k20*log(k20))/(2*k20)-0.3*m0;%求不变假设下的能量函数 U1=((log(u)-k11).^2-k21*log(k21))/(2*k21)-0.3*m1;%求变化假设下的能量函数 bianhua3=(U0>U1);%第一次得到的变化图像 chazhi=1; while chazhi>3e-3 p0=exp(-U0)./(exp(-U0)+exp(-U1)); p1=exp(-U1)./(exp(-U0)+exp(-U1)); w0=(p0>p1).*p0; w1=(p1>p0).*p1; k10=sum(sum(w0.*log(u)))/(sum(sum(w0)));%求不变部分的均值 k11=sum(sum(w1.*log(u)))/(sum(sum(w1)));%求变化部分的均值 k20=sum(sum(((log(u)-k10).^2).*w0))/(sum(sum(w0)));%求不变部分的方差 k21=sum(sum(((log(u)-k11).^2).*w1))/(sum(sum(w1)));%求变化部分的方差 disp('求步长'); buchang=newton(w0,w1,m0,m1);%牛顿-拉夫逊算法求步长因子 XX1=wextend('2D','sym',bianhua3,[1,1]); for i=2:hang+1 for j=2:lie+1 window3=XX1(i-1:i+1,j-1:j+1); if XX1(i,j)==0 m0(i-1,j-1)=8-sum(sum(window3)); m1(i-1,j-1)=sum(sum(window3)); else m0(i-1,j-1)=9-sum(sum(window3)); m1(i-1,j-1)=sum(sum(window3))-1; end end end U0=((log(u)-k10).^2-k20*log(k20))/(2*k20)-buchang*m0; U1=((log(u)-k11).^2-k21*log(k21))/(2*k21)-buchang*m1; bianhua4=(U0>U1); chazhi=sum(sum(abs(bianhua3-bianhua4)))/(hang*lie); fprintf('连续两次的结果差值=%i\n',chazhi) bianhua3=bianhua4; % figure; % imshow(bianhua1); end figure;imshow(bianhua3,[]);title('正变化分割结果'); % b = bwmorph(bianhua1,'clean',Inf);%去除孤立点b % x3 = im2double(b); % B1 = ones(7,7); % cs1 = conv2(x3,B1,'same'); % cs1 = x3.*cs1; % for c3 = 2 : hang % for d3 = 2:lie % if cs1(c3,d3) < 24 % cs1(c3,d3) = 0; % else % cs1(c3,d3) = 1; % end % end % end bianhua1=bianhua1|bianhua3; figure; imshow(bianhua1); title('变化图像11111111'); imwrite(bianhua1,'12.png'); toc