www.pudn.com > code.zip > mainfunction.m, change:2015-03-15,size:10646b


%%  
% ############################## 
%           2014.5 
%   初始化,读入图片,并进行滤波 
% ############################## 
clc; 
close all; 
tic; 
N=3; 
I=imread('1999.04.bmp'); 
Ic=imread('1999.05.bmp'); 
Ib=I(:,:,1);                 % 取出灰度图 
Icb=Ic(:,:,1);               % Ic表示和I进行比较(compare) 
Ibb=lee(Ib);               % 进行frost滤波  偶尔会出现out of memory的提示!(211_1.bmp图片大小是257*285) 
Icbb=lee(Icb);             % 相当于进行预处理 
I=double(Ibb);               % 为什么要加double?因为图片一般是uint8类型的! 
Ic=double(Icbb);             % 转为double类型后为下面的比值打下基础 
Dif=abs(log(I./Ic));         % 构造对数比值法构造差异图,作为reference map,是double类型 
figure(1); 
subplot(1,2,1),imshow(Ibb);  % 如果是imshow(I)则结果会出现一片空白,因为是double类型 
title('图1的lee滤波图');     
subplot(1,2,2),imshow(Icbb); 
title('图2的lee滤波图'); 
figure(2); 
imshow(Dif);                 % 这里加上uint8(Dif)则出现的图像是全黑的,因此这里不加uint8 
title('两幅图的对数比值差异图'); 
[M1,N1]=size(I); 
[M2,N2]=size(Ic); 
 
%% 
%  ############################### 
%     图像I的N个像素边界镜像拓展 
%  ############################### 
Id=zeros(M1+2*N,(N1+2*N)); 
for i=1:N 
    for j=1:N 
        Id(i,j)=I(N-i+1,N-j+1);  %将原图片I的左上角N*N的正方形块中心对称到Id的左上角上 
    end 
end 
 
for i=(M1+2*N-N+1):(M1+2*N) 
    for j=1:N 
        Id(i,j)=I(2*M1+N+1-i,N-j+1);  %将原图片I的左下角N*N的正方形块中心对称到Id的左下角上 
    end 
end 
 
for i=1:N 
    for j=(N1+2*N-N+1):(N1+2*N) 
        Id(i,j)=I(N-i+1,2*N1+N+1-j);  %将原图片I的右上角N*N的正方形块中心对称到Id的右上角上 
    end 
end 
 
for i=(M1+2*N-N+1):(M1+2*N) 
    for j=(N1+2*N-N+1):(N1+2*N) 
        Id(i,j)=I(2*M1+N+1-i,2*N1+N+1-j);  %将原图片I的右下角N*N的正方形块中心对称到Id的右下角上 
    end 
end 
 
for i=1:N 
    for j=(N+1):(N1+2*N-N) 
        Id(i,j)=I(N-i+1,j-N);  %将原图片I的上边部分N*N1的矩形块轴对称到Id的上边部分上 
    end 
end 
 
for i=(N+1):(M1+2*N-N) 
    for j=1:N 
        Id(i,j)=I(i-N,N-j+1);  %将原图片I的左边部分M1*N的矩形块轴对称到Id的左边部分上 
    end 
end 
 
for i=(N+1):(M1+2*N-N) 
    for j=(N1+2*N-N+1):(N1+2*N) 
        Id(i,j)=I(i-N,2*N1+N+1-j);  %将原图片I的右边部分M1*N的矩形块轴对称到Id的右边部分上 
    end 
end 
 
for i=(M1+2*N-N+1):(M1+2*N) 
    for j=(N+1):(N1+2*N-N) 
        Id(i,j)=I(2*M1+N+1-i,j-N);  %将原图片I的下边部分N*N1的矩形块轴对称到Id的下边部分上 
    end 
end 
 
for i=(N+1):(M1+2*N-N) 
    for j=(N+1):(N1+2*N-N) 
        Id(i,j)=I(i-N,j-N);    %将原图片I整体M1*N1复制到Id的中间,是得Id就是在I基础上加了一个“边框” 
    end 
end 
 
Id=uint8(Id);       %如果是imshow(Id)则会出现空白,需要加上uint8 
figure(3); 
imshow(Id);         %如果是imshow(Id)则会出现空白,需要加上uint8 
title('原图I的N像素拓展图像'); 
 
%% 
%  ############################### 
%     图像Ic的N个像素边界镜像拓展 
%  ############################### 
Icd=zeros(M2+2*N,(N2+2*N)); 
for i=1:N 
    for j=1:N 
        Icd(i,j)=Ic(N-i+1,N-j+1);  %将原图片Ic的左上角N*N的正方形块中心对称到Icd的左上角上 
    end 
end 
 
for i=(M2+2*N-N+1):(M2+2*N) 
    for j=1:N 
        Icd(i,j)=Ic(2*M2+N+1-i,N-j+1);  %将原图片Ic的左下角N*N的正方形块中心对称到Icd的左下角上 
    end 
end 
 
for i=1:N 
    for j=(N2+2*N-N+1):(N2+2*N) 
        Icd(i,j)=Ic(N-i+1,2*N2+N+1-j);  %将原图片Ic的右上角N*N的正方形块中心对称到Icd的右上角上 
    end 
end 
 
for i=(M2+2*N-N+1):(M2+2*N) 
    for j=(N2+2*N-N+1):(N2+2*N) 
        Icd(i,j)=Ic(2*M2+N+1-i,2*N2+N+1-j);  %将原图片Ic的右下角N*N的正方形块中心对称到Icd的右下角上 
    end 
end 
 
for i=1:N 
    for j=(N+1):(N2+2*N-N) 
        Icd(i,j)=Ic(N-i+1,j-N);  %将原图片Ic的上边部分N*N2的矩形块轴对称到Icd的上边部分上 
    end 
end 
 
for i=(N+1):(M2+2*N-N) 
    for j=1:N 
        Icd(i,j)=Ic(i-N,N-j+1);  %将原图片Ic的左边部分M2*N的矩形块轴对称到Icd的左边部分上 
    end 
end 
 
for i=(N+1):(M2+2*N-N) 
    for j=(N2+2*N-N+1):(N2+2*N) 
        Icd(i,j)=Ic(i-N,2*N2+N+1-j);  %将原图片Ic的右边部分M2*N的矩形块轴对称到Icd的右边部分上 
    end 
end 
 
for i=(M2+2*N-N+1):(M2+2*N) 
    for j=(N+1):(N2+2*N-N) 
        Icd(i,j)=Ic(2*M1+N+1-i,j-N);  %将原图片Ic的下边部分N*N2的矩形块轴对称到Icd的下边部分上 
    end 
end 
 
for i=(N+1):(M2+2*N-N) 
    for j=(N+1):(N2+2*N-N) 
        Icd(i,j)=Ic(i-N,j-N);    %将原图片Ic整体M2*N2复制到Icd的中间,使得Icd就是在Ic基础上加了一个“边框” 
    end 
end 
 
Icd=uint8(Icd);           %如果是imshow(Icd)则会出现空白,需要加上uint8 
figure(4); 
imshow(Icd);              %如果是imshow(Icd)则会出现空白,需要加上uint8 
title('原图Ic的N像素拓展图像'); 
 
%%  
%  ############################# 
%     构造像素(i,j)的特征向量 
%  ############################# 
I1=Id(N:(M1+2*N-N+1),N:(N1+2*N-N+1));  %由于是取N*N的矩阵,而该函数是从最开始的行列算起,故去掉最外围的(N-1)行(N-1)列 
II=I1';                                %由于im2col函数的特殊性,需要转置 
III=zeros(N^2,1,M1*N1); 
III1=im2col(II,[N,N]);                 %一共得到M1*N1个列向量,也就是yij向量(每个像素的特征向量) 
for p=1:M1*N1 
    III(:,:,p)=III1(:,p); 
end 
%disp(III); 
%y1ij=III(:,i);                   %取出矩阵第i列的列向量 
 
Ic1=Icd(N:(M2+2*N-N+1),N:(N2+2*N-N+1));    %由于是取N*N的矩阵,而该函数是从最开始的行列算起,故去掉最外围的两行两列 
IIc=Ic1';                         %由于im2col函数的特殊性,需要转置 
IIIc=zeros(N^2,1,M2*N2); 
IIIc1=im2col(IIc,[N,N]);           %一共得到M2*N2个列向量,也就是yij向量(每个像素的特征向量) 
for p=1:M2*N2 
    IIIc(:,:,p)=IIIc1(:,p); 
end 
%disp(IIIc); 
%toc; 
%y1ij=III(:,i);                 %取出矩阵第i列的列向量 
 
%%  
%  #################################### 
%     构造两幅图像像素(i,j)的局部字典 
%  ####################################  
D1=zeros(N^2,(N+2)^2,M1*N1);             %共有M1*N1个字典D 
Id1=zeros((2*N+1),(2*N+1),M1*N1);         %Id1是过渡矩阵 
for i=1:M1                     %由于是(2*N+1)*(2*N+1)的搜索窗,故中心像素应该是从除去外围的三行三列开始的 
    for j=1:N1 
        Id1(:,:,((i-1)*N1+j))=Id(i:(i+2*N),j:(j+2*N)); 
        Id1(:,:,((i-1)*N1+j))=Id1(:,:,((i-1)*(N1-2*N)+j))';             %转置 
        D1(:,:,((i-1)*N1+j))=im2col(Id1(:,:,((i-1)*N1+j)),[N,N]);   % ((i-1)*(M1-2*N)+j)表示像素(i,j)对应第((i-1)*(M1-2*N)+j)个字典         
    end 
end 
 
D2=zeros(N^2,(N+2)^2,M2*N2);             %共有M2*N2个字典D 
Id2=zeros((2*N+1),(2*N+1),M2*N2); 
for i=1:M2               %由于是(2*N+1)*(2*N+1)的搜索窗,故中心像素应该是从除去外围的三行三列开始的 
    for j=1:N2 
        Id2(:,:,((i-1)*N2+j))=Icd(i:(i+2*N),j:(j+2*N)); 
        Id2(:,:,((i-1)*N2+j))=Id2(:,:,((i-1)*N2+j))'; 
        D2(:,:,((i-1)*N2+j))=im2col(Id2(:,:,((i-1)*N2+j)),[N,N]);   % ((i-1)*(M1-2*N)+j)表示像素(i,j)对应第((i-1)*(M1-2*N)+j)个字典 
    end 
end 
 
%% 
% ############################## 
%   利用OMP算法求出稀疏表示系数 
% ############################## 
Q=N^2; 
A=zeros((N+2)^2,1,M1*N1); 
Ac=zeros((N+2)^2,1,M2*N2); 
III=double(III); 
IIIc=double(IIIc); 
D1=double(D1); 
% for i=1:(M1-2*N)*(N1-2*N)                      %将字典D1、D2的原子进行归一化 
%    for k=1:(N+2)^2 
%         D1(:,k,i)=D1(:,k,i)/norm(D1(:,k,i)); 
%     end 
% end 
D2=double(D2); 
% for i=1:(M1-2*N)*(N1-2*N) 
%     for k=1:(N+2)^2 
%        D2(:,k,i)=D2(:,k,i)/norm(D2(:,k,i)); 
%     end 
% end 
 
for i=1:M1*N1 
    [A(:,:,i)]=OMP(D2(:,:,i),III(:,:,i),Q);     
end 
for k=1:M2*N2 
    [Ac(:,:,k)]=OMP(D1(:,:,k),IIIc(:,:,k),Q);  
end 
 
%% 
% ######################## 
%   求解L1范并构造差异图 
% ######################## 
L1=zeros(M1,N1);        
L1c=zeros(M2,N2);  
% ===================================================== 
%         rem(n,m)或mod(n,m)表示n/m的余数 
%         圆整和求余函数(Rounding and remainder) 
%         ceil 朝正无穷大方向取整 
%         fix 朝零方向取整 
%         floor 朝负无穷大方向取整 
%         mod 模数求余  用法是mod(a,b) 
%         rem 求余数 
%         round 四舍五入取整 
%         sign 符号函数 
%====================================================== 
for i=1:M1*N1 
    p=ceil(i/N1); 
    if mod(i,N1)==0 
        q=fix(i/N1)*N1; 
    else 
        q=mod(i,N1); 
    end 
    L1(p,q)=norm(A(:,:,i),1);    %n = norm(X,1)    %求1-范数 
end 
for k=1:M2*N2 
    r=ceil(k/N1); 
    if mod(k,N1)==0 
        s=fix(k/N1)*N1; 
    else 
        s=mod(k,N1); 
    end 
    L1c(r,s)=norm(Ac(:,:,k),1);    
end 
 
L=zeros(M1,N1);  
for i=1:M1 
    for j=1:N1 
        L(i,j)=abs(L1(i,j)-L1c(i,j));     %构造差异图像        
    end 
end 
L=double(L); 
 figure; 
 imshow(L); 
 title('两幅图像的差异图像'); 
 
%% 
% ################################### 
%      对差异图像进行最大阈值分割 
% ################################### 
LL=L; 
count=imhist(LL); 
[m,n]=size(LL); 
S=m*n; 
K=256; 
count=count/S;   %%每一个像素的分布概率 
count; 
 
for i=1:K 
    if count(i)~=0 
        st=i-1; 
        break; 
    end 
end 
st; 
for i=K:-1:1 
    if count(i)~=0 
        nd=i-1; 
        break; 
    end 
end 
nd; 
f=count(st+1:nd+1);  %f是每个灰度出现的概率 
size(f) 
E=[]; 
for Th=st:nd-1   %%%设定初始分割阈值为Th 
    av1=0; 
    av2=0; 
    Pth=sum(count(1:Th+1)); 
%%%第一类的平均相对熵为 
    for i=0:Th 
        av1=av1-count(i+1)/Pth*log(count(i+1)/Pth+0.00001); 
    end 
%%%第二类的平均相对熵为 
    for i=Th+1:K-1 
        av2=av2-count(i+1)/(1-Pth)*log(count(i+1)/(1-Pth)+0.00001); 
    end 
       E(Th-st+1)=av1+av2; 
end 
position=find(E==(max(E)));  %找出最大值的位置 
th=st+position-1; 
 
 for i=1:m 
    for j=1:n 
        if LL(i,j)>th 
            LL(i,j)=255; 
        else 
            LL(i,j)=0; 
        end 
    end 
end  
figure; 
imshow(LL); 
toc; 
 
%%   
% ######################################### 
%       统计变化点像素和非变化点像素个数 
% ######################################### 
index1=find(abs(L)>10); 
b=size(index1);               %  取出变化像素点的个数 
index2=find(abs(L)<0.1); 
c=size(index2);               %  取出像素点不变化的点的个数 
 
fprintf('变化点像素个数是%d\n',b(1,1)); 
fprintf('非变化点像素个数是%d\n',c(1,1)); 
 
% ################ 
%    统计错误率 
% ################ 
ref=imread('reference_Bern.bmp'); 
ref=ref(:,:,1);                 %这一句不可缺省,否则会出现“??? Array dimensions must match for binary arrayop.”的提示 
I11=uint8(L);                   %加了FCM算法后,错误率很高达38%了! 
I22=ref; 
figure; 
subplot(1,2,1),imshow(I22); 
title('两幅图像的变化参考图'); 
subplot(1,2,2),imshow(L);       %如果是显示I11,则是全黑,所以不加uint8 
title('两幅图像的差异图像'); 
 
[N,LLL]=size(I); 
[T,J]=cuowu(I11,I22,N,LLL);     %统计误检及漏检点数 
K=T+J; 
disp(sprintf('结果漏检%d个像素点,误检%d个像素点',T,J)); 
disp(sprintf('结果图像共有%d个像素,其中参考误差像素个数为%d',N*LLL,K));  %也可以用fprintf命令 
K=K/(N*LLL)*100; 
disp(sprintf('此方法与参考图像的误差率为%2.4f个百分点',K));