ms.m

图像抽取及误差控制

function  [kkkk,fff,ff4,f,emse,PSNR]=ms
%ERROR CONCTROL Summary of this function goes here
%   Detailed explanation goes here
%读取图像
f=imread('LENA.BMP');

f=single(f);
%第一次抽取图像
for i=1:256
for j=1:256
f3(i,j)=f(2*i,2*j);
end
end;
%第二次抽取图像
for i=1:128
for j=1:128
f2(i,j)=f3(2*i,2*j);
end
end;
%第三次抽取图像
for i=1:64
for j=1:64
f1(i,j)=f2(2*i,2*j);
end
end;
%第四次抽取图像
for i=1:32
for j=1:32
f0(i,j)=f1(2*i,2*j);
end
end;
%预测插值
for i=3:255
for j=3:255
da1(i,j)=f(2*i-1,2*j)-9/16*(f3(i,j)+f3(i-1,j))+1/16*(f3(i+1,j)+f3(i-2,j));
da2(i,j)=f(2*i,2*j-1)-9/16*(f3(i,j)+f3(i,j-1))+1/16*(f3(i,j+1)+f3(i,j-2));
da3(i,j)=f(2*i-1,2*j-1)-(9/16*(9/16*(f3(i,j)+f3(i-1,j)+f3(i,j-1)+f3(i-1,j-1))...
-1/16*(f3(i,j+1)+f3(i-1,j+1)+f3(i,j-2)+f3(i-1,j-2)))-1/16*(9/16*(f3(i+1,j)...
+f3(i-2,j)+f3(i+1,j-1)+f3(i-2,j-1))-1/16*(f3(i+1,j+1)+f3(i-2,j+1)+f3(i+1,j-2)...
+f3(i-2,j-2))));
end
end;

for i=3:127
for j=3:127
db1(i,j)=f3(2*i-1,2*j)-9/16*(f2(i,j)+f2(i-1,j))+1/16*(f2(i+1,j)+f2(i-2,j));
db2(i,j)=f3(2*i,2*j-1)-9/16*(f2(i,j)+f2(i,j-1))+1/16*(f2(i,j+1)+f2(i,j-2));
db3(i,j)=f3(2*i-1,2*j-1)-(9/16*(9/16*(f2(i,j)+f2(i-1,j)+f2(i,j-1)+f2(i-1,j-1))...
-1/16*(f2(i,j+1)+f2(i-1,j+1)+f2(i,j-2)+f2(i-1,j-2)))-1/16*(9/16*(f2(i+1,j)...
+f2(i-2,j)+f2(i+1,j-1)+f2(i-2,j-1))-1/16*(f2(i+1,j+1)+f2(i-2,j+1)...
+f2(i+1,j-2)+f2(i-2,j-2))));
end
end;

for i=3:63
for j=3:63
dc1(i,j)=f2(2*i-1,2*j)-9/16*(f1(i,j)+f1(i-1,j))+1/16*(f1(i+1,j)+f1(i-2,j));
dc2(i,j)=f2(2*i,2*j-1)-9/16*(f1(i,j)+f1(i,j-1))+1/16*(f1(i,j+1)+f1(i,j-2));
dc3(i,j)=f2(2*i-1,2*j-1)-(9/16*(9/16*(f1(i,j)+f1(i-1,j)+f1(i,j-1)+f1(i-1,j-1))...
-1/16*(f1(i,j+1)+f1(i-1,j+1)+f1(i,j-2)+f1(i-1,j-2)))-1/16*(9/16*(f1(i+1,j)...
+f1(i-2,j)+f1(i+1,j-1)+f1(i-2,j-1))-1/16*(f1(i+1,j+1)+f1(i-2,j+1)+f1(i+1,j-2)...
+f1(i-2,j-2))));
end
end;

for i=3:31
for j=3:31
dd1(i,j)=f1(2*i-1,2*j)-9/16*(f0(i,j)+f0(i-1,j))+1/16*(f0(i+1,j)+f0(i-2,j));
dd2(i,j)=f1(2*i,2*j-1)-9/16*(f0(i,j)+f0(i,j-1))+1/16*(f0(i,j+1)+f0(i,j-2));
dd3(i,j)=f1(2*i-1,2*j-1)-(9/16*(9/16*(f0(i,j)+f0(i-1,j)+f0(i,j-1)+f0(i-1,j-1))...
-1/16*(f0(i,j+1)+f0(i-1,j+1)+f0(i,j-2)+f0(i-1,j-2)))-1/16*(9/16*(f0(i+1,j)...
+f0(i-2,j)+f0(i+1,j-1)+f0(i-2,j-1))-1/16*(f0(i+1,j+1)+f0(i-2,j+1)...
+f0(i+1,j-2)+f0(i-2,j-2))));
end
end;

%差值转换为整数
da11=round(da1);
da12=round(da2);
da13=round(da3);
db11=round(db1);
db12=round(db2);
db13=round(db3);
dc11=round(dc1);
dc12=round(dc2);
dc13=round(dc3);
dd11=round(dd1);
dd12=round(dd2);
dd13=round(dd3);
%再把插值转换为单精度类型
da11=single(da11);
da12=single(da12);
da13=single(da13);
db11=single(db11);
db12=single(db12);
db13=single(db13);
dc11=single(dc11);
dc12=single(dc12);
dc13=single(dc13);
dd11=single(dd11);
dd12=single(dd12);
dd13=single(dd13);
%for i=3:255
%for j=3:255
%da11(i,j)=2*e*round(da1(i,j)/2/e);
%da12(i,j)=2*e*round(da2(i,j)/2/e);
%da13(i,j)=2*e*round(da3(i,j)/2/e);
%end
%end

%for i=3:127
%for j=3:127
%db11(i,j)=2*e*round(db1(i,j)/2/e);
%db12(i,j)=2*e*round(db2(i,j)/2/e);
%db13(i,j)=2*e*round(db3(i,j)/2/e);
%end
%end

%for i=3:63
%for j=3:63
%dc11(i,j)=2*e*round(dc1(i,j)/2/e);
%dc12(i,j)=2*e*round(dc2(i,j)/2/e);
%dc13(i,j)=2*e*round(dc3(i,j)/2/e);
%end
%end

%for i=3:31
%for j=3:31
%dd11(i,j)=2*e*round(dd1(i,j)/2/e);
%dd12(i,j)=2*e*round(dd2(i,j)/2/e);
%dd13(i,j)=2*e*round(dd3(i,j)/2/e);
%end
%end




%图像的重构
for i=1:32
for j=1:32
ff1(2*i,2*j)=f0(i,j);
end
end
for i=1:4
for j=1:64
ff1(i,j)=f1(i,j);
end
end;
for j=1:4
for i=1:64
ff1(i,j)=f1(i,j);
end
end;
for j=1:64
ff1(63,j)=f1(63,j);
ff1(64,j)=f1(64,j);
end
for i=1:64
ff1(i,63)=f1(i,63);
ff1(i,64)=f1(i,64);
end
for i=3:31
for j=3:31
ff1(2*i-1,2*j)=9/16*(f0(i,j)+f0(i-1,j))-1/16*(f0(i+1,j)+f0(i-2,j))+dd11(i,j);
ff1(2*i,2*j-1)=9/16*(f0(i,j)+f0(i,j-1))-1/16*(f0(i,j+1)+f0(i,j-2))+dd12(i,j);
ff1(2*i-1,2*j-1)=(9/16*(9/16*(f0(i,j)+f0(i-1,j)+f0(i,j-1)+f0(i-1,j-1))-1/16*(f0(i,j+1)+f0(i-1,j+1)+f0(i,j-2)+f0(i-1,j-2)))-1/16*(9/16*(f0(i+1,j)+f0(i-2,j)+f0(i+1,j-1)+f0(i-2,j-1))-1/16*(f0(i+1,j+1)+f0(i-2,j+1)+f0(i+1,j-2)+f0(i-2,j-2))))+dd13(i,j);
end
end




for i=1:64
for j=1:64
ff2(2*i,2*j)=ff1(i,j);
end
end
for i=1:4
for j=1:128
ff2(i,j)=f2(i,j);
end
end;
for j=1:4
for i=1:128
ff2(i,j)=f2(i,j);
end
end;
for j=1:128
ff2(127,j)=f2(127,j);
ff2(128,j)=f2(128,j);
end
for i=1:128
ff2(i,127)=f2(i,127);
ff2(i,128)=f2(i,128);
end
for i=3:63
for j=3:63
ff2(2*i-1,2*j)=9/16*(ff1(i,j)+ff1(i-1,j))-1/16*(ff1(i+1,j)+ff1(i-2,j))+dc11(i,j);
ff2(2*i,2*j-1)=9/16*(ff1(i,j)+ff1(i,j-1))-1/16*(ff1(i,j+1)+ff1(i,j-2))+dc12(i,j);
ff2(2*i-1,2*j-1)=(9/16*(9/16*(ff1(i,j)+ff1(i-1,j)+ff1(i,j-1)+ff1(i-1,j-1))...
-1/16*(ff1(i,j+1)+ff1(i-1,j+1)+ff1(i,j-2)+ff1(i-1,j-2)))-1/16*(9/16*(ff1(i+1,j)...
+ff1(i-2,j)+ff1(i+1,j-1)+ff1(i-2,j-1))-1/16*(ff1(i+1,j+1)+ff1(i-2,j+1)...
+ff1(i+1,j-2)+ff1(i-2,j-2))))+dc13(i,j);
end
end




for i=1:128
for j=1:128
ff3(2*i,2*j)=ff2(i,j);
end
end
for i=1:4
for j=1:256
ff3(i,j)=f3(i,j);
end
end;
for j=1:4
for i=1:256
ff3(i,j)=f3(i,j);
end
end;
for j=1:256
ff3(255,j)=f3(255,j);
ff3(256,j)=f3(256,j);
end
for i=1:256
ff3(i,255)=f3(i,255);
ff3(i,256)=f3(i,256);
end
for i=3:127
for j=3:127
ff3(2*i-1,2*j)=9/16*(ff2(i,j)+ff2(i-1,j))-1/16*(ff2(i+1,j)+ff2(i-2,j))+db11(i,j);
ff3(2*i,2*j-1)=9/16*(ff2(i,j)+ff2(i,j-1))-1/16*(ff2(i,j+1)+ff2(i,j-2))+db12(i,j);
ff3(2*i-1,2*j-1)=(9/16*(9/16*(ff2(i,j)+ff2(i-1,j)+ff2(i,j-1)+ff2(i-1,j-1))...
-1/16*(ff2(i,j+1)+ff2(i-1,j+1)+ff2(i,j-2)+ff2(i-1,j-2)))-1/16*(9/16*...
(ff2(i+1,j)+ff2(i-2,j)+ff2(i+1,j-1)+ff2(i-2,j-1))-1/16*(ff2(i+1,j+1)...
+ff2(i-2,j+1)+ff2(i+1,j-2)+ff2(i-2,j-2))))+db13(i,j);
end
end


for i=1:256
for j=1:256
ff4(2*i,2*j)=ff3(i,j);
end
end
for i=1:4
for j=1:512
ff4(i,j)=f(i,j);
end
end;
for j=1:4
for i=1:512
ff4(i,j)=f(i,j);
end
end;
for j=1:512
ff4(511,j)=f(511,j);
ff4(512,j)=f(512,j);
end
for i=1:512
ff4(i,511)=f(i,511);
ff4(i,512)=f(i,512);
end
for i=3:255
for j=3:255
ff4(2*i-1,2*j)=9/16*(ff3(i,j)+ff3(i-1,j))-1/16*(ff3(i+1,j)+ff3(i-2,j))+da11(i,j);
ff4(2*i,2*j-1)=9/16*(ff3(i,j)+ff3(i,j-1))-1/16*(ff3(i,j+1)+ff3(i,j-2))+da12(i,j);
ff4(2*i-1,2*j-1)=(9/16*(9/16*(ff3(i,j)+ff3(i-1,j)+ff3(i,j-1)+ff3(i-1,j-1))-...
1/16*(ff3(i,j+1)+ff3(i-1,j+1)+ff3(i,j-2)+ff3(i-1,j-2)))-1/16*(9/16*(ff3(i+1,j)...
+ff3(i-2,j)+ff3(i+1,j-1)+ff3(i-2,j-1))-1/16*(ff3(i+1,j+1)+ff3(i-2,j+1)...
+ff3(i+1,j-2)+ff3(i-2,j-2))))+da13(i,j);
end
end
%ff4为重构图像
%fff为重构图像与原图像的差矩阵
fig=uint8(f);
ff4=uint8(ff4);
fff=fig-ff4;
figure,imshow(fig)
figure,imshow(ff4)

%误差分析

ffff=reshape(fff,1,512*512);
ffff=single(ffff);
emse=(var(ffff,1));
PSNR=10*(log(512*512*256*256/sum(ffff.*ffff)));










for i=1:253
for j=1:253
daf1(i,j)=da11(i+2,j+2);
daf2(i,j)=da12(i+2,j+2);
daf3(i,j)=da13(i+2,j+2);
end
end;

for i=1:125
for j=1:125
dbf1(i,j)=db11(i+2,j+2);
dbf2(i,j)=db12(i+2,j+2);
dbf3(i,j)=db13(i+2,j+2);
end
end;

for i=1:61
for j=1:61
dcf1(i,j)=dc11(i+2,j+2);
dcf2(i,j)=dc12(i+2,j+2);
dcf3(i,j)=dc13(i+2,j+2);
end
end;

for i=1:29
for j=1:29
ddf1(i,j)=dd11(i+2,j+2);
ddf2(i,j)=dd12(i+2,j+2);
ddf3(i,j)=dd13(i+2,j+2);
end
end;





a=reshape(daf1,1,numel(daf1));
b=reshape(daf2,1,numel(daf2));
c=reshape(daf3,1,numel(daf3));
d=reshape(dbf1,1,numel(dbf1));
e=reshape(dbf2,1,numel(dbf2));
fg=reshape(dbf3,1,numel(dbf3));
g=reshape(dcf2,1,numel(dcf2));
h=reshape(dcf1,1,numel(dcf1));
ii=reshape(dcf3,1,numel(dcf3));
k=reshape(ddf1,1,numel(ddf1));
l=reshape(ddf2,1,numel(ddf2));
m=reshape(ddf3,1,numel(ddf3));
n=reshape(f0,1,numel(f0));
qq0=[a b c d e fg g h ii  k l m n];
% 
qq1=[f1(1,:) f1(3,:) f1(63,:)];



for i=1:32
qq2(1,i)=f1(2,2*i);
end;

for i=1:32
qq3(1,i)=f1(4,2*i);
end;

for i=1:32
qq4(1,i)=f1(64,2*i);
end;





for i=1:58
qq5(1,i)=f1(i+4,1);
qq6(1,i)=f1(i+4,3);
qq7(1,i)=f1(i+4,63);
end;


for i=1:29
qq8(1,i)=f1(2*(i+1)+1,2);
qq9(1,i)=f1(2*(i+1)+1,4);
qq10(1,i)=f1(2*(i+1)+1,64);

end

qqq=[qq0 qq1 qq2 qq3 qq4 qq5 qq6 qq7 qq8 qq9 qq10];




qq11=[f2(1,:) f2(3,:) f2(127,:)];

for i=1:64
qq12(1,i)=f2(2,2*i);
end;

for i=1:64
qq13(1,i)=f2(4,2*i);
end;

for i=1:64
qq14(1,i)=f2(128,2*i);
end;





for i=1:122
qq15(1,i)=f2(i+4,1);
qq16(1,i)=f2(i+4,3);
qq17(1,i)=f2(i+4,63);
end;


for i=1:61
qq18(1,i)=f2(2*(i+1)+1,2);
qq19(1,i)=f2(2*(i+1)+1,4);
qq110(1,i)=f2(2*(i+1)+1,128);

end

qqq1=[qq11 qq12 qq13 qq14 qq15 qq16 qq17 qq18 qq19 qq110];





qq21=[f3(1,:) f3(3,:) f3(255,:)];

for i=1:128
qq22(1,i)=f3(2,2*i);
end;

for i=1:128
qq23(1,i)=f3(4,2*i);
end;

for i=1:128
qq24(1,i)=f3(256,2*i);
end;





for i=1:250
qq25(1,i)=f3(i+4,1);
qq26(1,i)=f3(i+4,3);
qq27(1,i)=f3(i+4,255);
end;


for i=1:125
qq28(1,i)=f3(2*(i+1)+1,2);
qq29(1,i)=f3(2*(i+1)+1,4);
qq210(1,i)=f3(2*(i+1)+1,256);

end

qqq2=[qq21 qq22 qq23 qq24 qq25 qq26 qq27 qq28 qq29 qq210];
%  ggg=[qqq qqq1 qqq2 ];
qq31=[f(1,:) f(3,:) f(511,:)];

for i=1:256
qq32(1,i)=f(2,2*i);
end;

for i=1:256
qq33(1,i)=f(4,2*i);
end;

for i=1:256
qq34(1,i)=f(512,2*i);
end;





for i=1:506
qq35(1,i)=f(i+4,1);
qq36(1,i)=f(i+4,3);
qq37(1,i)=f(i+4,255);
end;


for i=1:253
qq38(1,i)=f(2*(i+1)+1,2);
qq39(1,i)=f(2*(i+1)+1,4);
qq310(1,i)=f(2*(i+1)+1,512);

end

qqq3=[qq31 qq32 qq33 qq34 qq35 qq36 qq37 qq38 qq39 qq310];

kkk=[qqq qqq1 qqq2 qqq3];
kkkk=reshape(kkk,512,512);