1 内容介绍
网络技术和多媒体技术迅猛发展,为了更好地保障图像信息传输的安全性和可靠性,解决数字图像的版权保护问题,图像信息隐藏技术已成为图像处理领域的研究热点之一.采用MATLAB开发环境,提出了一种图像信息隐藏算法,实现了图像信息隐藏的两大基本功能,即信息的嵌入和信息的提取.同时,模拟常见的攻击方法,验证了添加隐藏信息后的图像的抗攻击能力.实验表明,该算法具有较强的鲁棒性,添加隐藏信息后的图像能有效地应对多种常见攻击,且隐藏的信息具有良好的不可见性.
2 仿真代码
I2_2=attack1(d,I0_wm,fig);
%%%去除权威0的行和列，使本文算法可以抵抗RT、RST的组合攻击
% [ro2,co2]=size(I2_2);
% m=0;I2_21=[];
% for i=1:ro2
% if sum(uint8(I2_2(i,:)))>512
% m=m+1;
% I2_21(m,:)=uint8(I2_2(i,:));
% end
% end
% n=0;
% for i=1:co2
% if sum(I2_21(:,i))>512
% n=n+1;
% I2_22(:,n)=I2_21(:,i);
% end
% end
% I2_2=I2_22;
% figure,imshow(uint8(I2_2));title('去除行列为0后的图像')
% I1=im2double(imread('TestImages/lena.bmp'));
% rect=[20 20 460 460];
% I2_2=imcrop(uint8(I2_0),rect);
% 先scale后rotation
% Scl=614;
% I2_1=imresize(I2_0,[Scl Scl],'bilinear'); %%无论是imrotate还是imresize采用'bicubic'是最好的
% I2_2=imrotate(I2_1,10,'bilinear');
I2=im2double(uint8(I2_2)); %先必须化为uint8型，然后用im2double归一化
% Get the Key Points
%% 受攻击后的图像与原图像（备注：此处有误，其实应该为嵌入水印的载体图像）进行第一次特征点的匹配
Options.upright=true;
Options.tresh=0.0002;
Ipts1=OpenSurf(I1,Options);
Ipts2=OpenSurf(I2,Options);
% 挑选Scale在一定范围内的点，因尺度太大或太小，受攻击后易丢失，但受攻击后的图像不需要选择
% m=0;
% for i=1:length(Ipts1)
% if Ipts1(i).scale>=2.0&&Ipts1(i).scale<=6.0
% m=m+1;
% Ipts1_impv(m)=Ipts1(i);
% end
% end
% save('fp_mult.mat','Ipts1_impv');
% fp_tst=load('fp_mult.mat');
% Ipts1_impv=fp_tst.Ipts1_impv;
point_cnt=0;ip_tmp1=[];
for i=1:length(Ipts1)
ip=Ipts1(i);
S = 2 * fix(2.5 * ip.scale);
% R = fix(S / 2);
if S>=10.0&&S<=30
ip_tmp1=[ip_tmp1;Ipts1(i)];
point_cnt=point_cnt+1;
% pt = [(ip.x), (ip.y)];
% ptR = [(R * cos(ip.orientation)), (R * sin(ip.orientation))];
%
% if(ip.laplacian >0), myPen =[0 0 1]; else myPen =[1 0 0]; end
%
% rectangle('Curvature', [1 1],'Position', [pt(1)-R, pt(2)-R, S, S],'EdgeColor',myPen);
%
% plot([pt(1), pt(1)+ptR(1)]+1,[pt(2), pt(2)+ptR(2)]+1,'g');
end
end
% orientation_max=max(ip_tmp1.orientation);
% nc_max=max(Nc_sta(:,1));
% tmp=Nc_sta(:,1)>=0.86*nc_max;
% nc_coord_alt=Nc_sta(tmp,:);
[a1,idx1]=sort(abs([ip_tmp1.doh]),'descend');
ip_sort=ip_tmp1(idx1);
nn=0;ip_fp=ip_sort(1);
for i=2:length(ip_sort)
if nn>=17
break;
end
S_sort = 2 * fix(2.5 * ip_sort(i).scale);
R_sort = fix(S_sort/2);
% mm=0;
kk=0;
for j=1:length(ip_fp)
S_fp = 2 * fix(2.5 * ip_fp(j).scale);
R_fp = fix(S_fp/2);
% mm=mm+1;
if sqrt((ip_fp(j).x-ip_sort(i).x)^2+(ip_fp(j).y-ip_sort(i).y)^2)>=R_sort+R_fp
kk=kk+1;
end
end
if kk==length(ip_fp)
nn=nn+1;
ip_fp=[ip_fp;ip_sort(i)];
end
end
save('fp_invT_480_lna1.mat','ip_fp');
figure,imshow(I1,[]);hold on;
for i=1:length(ip_fp)
ip=ip_fp(i);
S = 2 * fix(2.5 * ip.scale);
R = fix(S / 2);
% if S>=12&&S<=30
% ip_tmp1=[ip_tmp1;ipts(i)];
% point_cnt=point_cnt+1;
pt = [(ip.x), (ip.y)];
ptR = [(R * cos(ip.orientation)), (R * sin(ip.orientation))];
if(ip.laplacian >0), myPen =[0 0 1]; else myPen =[1 0 0]; end
rectangle('Curvature', [1 1],'Position', [pt(1)-R, pt(2)-R, S, S],'EdgeColor',myPen);
plot([pt(1), pt(1)+ptR(1)]+1,[pt(2), pt(2)+ptR(2)]+1,'g');
end
fp1=load('fp_invT_480_lna1.mat');
Ipts1_impv=fp1.ip_fp;
% Ipts2_impv=[];
n=0;
for i=1:length(Ipts2)
if Ipts2(i).scale>=2.0&&Ipts2(i).scale<=6.5
n=n+1;
Ipts2_impv(n)=Ipts2(i);
end
end
% Ipts2_impv=Ipts2;
% Put the landmark descriptors in a matrix
D1 = reshape([Ipts1_impv.descriptor],64,[]);
D2 = reshape([Ipts2_impv.descriptor],64,[]);
% Find the best
% matches，最近邻次近邻比值的阈值设定法挑选，Ransac算法再进行一次挑选，是否要用到最小二乘法拟合那个矩阵？？？
err=zeros(1,length(Ipts1_impv));
% cor1=1:length(Ipts1_impv);
cor2=zeros(1,length(Ipts1_impv));
ratio=zeros(length(Ipts1_impv),1);
dis_store=zeros(length(Ipts1_impv),2);
dis_ratio=zeros(length(Ipts1_impv),5);
for i=1:length(Ipts1_impv)
distance=sum((D2-repmat(D1(:,i),[1 length(Ipts2_impv)])).^2,1);
dis_sort=sort(distance,'ascend');
dis_store(i,:)=[dis_sort(1),dis_sort(2),];
ratio(i)=dis_sort(1)/dis_sort(2);
[err(i),cor2(i)]=min(distance);
dis_ratio(i,:)=[dis_store(i,:) ratio(i) i cor2(i)];
end
dis_ratio1=sortrows(dis_ratio,3);
mask1=dis_ratio1(:,3)<=0.9;
dis_ratio2=zeros(sum(mask1),5);
nn=0;
for i=1:length(dis_ratio1)
if mask1(i)
nn=nn+1;
dis_ratio2(nn,:)=dis_ratio1(i,:);
end
end
% 利用Ransac算法对匹配的点再进行一次挑选
max_itera=50; %设置最大迭代次数
sigma=1.8; %设置拟合矩阵[m11,m12,m13;m21,m22,m23;0,0,1]还原的偏差
pretotal=0; %符合拟合模型的数据的个数
k=0;
% sample=round(18*rand(50,3));
% save('Index.mat','sample');
sam_tmp=load('Index.mat');
Samp=sam_tmp.sample;
while pretotal <= size(dis_ratio2,1)*0.75 && k<max_itera %有2/3的数据符合拟合模型或达到最大迭代次数就可以退出了
SampIndex=Samp(k+1,:);
% SampIndex=round(1+(size(dis_ratio2,1)-1)*rand(3,1)); %产生三个随机索引，找样本用，floor向下取整
samp1_tmp=dis_ratio2(SampIndex,4); %对原数据随机抽样两个样本
samp2_tmp=dis_ratio2(SampIndex,5);
samp1=[Ipts1_impv(samp1_tmp).y;Ipts1_impv(samp1_tmp).x]; % 此处特别注意，其实在利用SURF算法时，其中有一步
samp2=[Ipts2_impv(samp2_tmp).y;Ipts2_impv(samp2_tmp).x]; % 通过旋转Haar小波模板求响应值，将X，Y的坐标对调了
att_matr=Matrix_solv(samp1,samp2); %对两组数据拟合出矩阵，或其他变种拟合方法
total=0;dis_ratio3=[];
for j=1:size(dis_ratio2,1)
x_ini=Ipts1_impv(dis_ratio2(j,4)).y;y_ini=Ipts1_impv(dis_ratio2(j,4)).x;
x_atk=Ipts2_impv(dis_ratio2(j,5)).y;y_atk=Ipts2_impv(dis_ratio2(j,5)).x;
tmp=att_matr*[x_atk;y_atk;1];
x_red=tmp(1);y_red=tmp(2);
if sqrt((x_ini-x_red)^2+(y_ini-y_red)^2)<sigma
total=total+1;
dis_ratio3=[dis_ratio3;dis_ratio2(j,:)];
end
end
% mask=abs(line*[data ones(size(data,1),1)]'); %求每个数据到拟合直线的距离
% total=sum(mask<sigma); %计算数据距离直线小于一定阈值的数据的个数
if total>pretotal %找到符合拟合矩阵数据最多的拟合矩阵
pretotal=total;
bestmatr=att_matr; %找到最好的拟合矩阵
best_disratio=dis_ratio3; %找到最符合条件的那些坐标
best_matr_r=att_matr;
end
k=k+1;
end
best_matr_r1=inv(best_matr_r);
% Sort matches on vector distance
% [err, ind]=sort(err);
cor1_1=best_disratio(:,4);
cor2_1=best_disratio(:,5);
%% 法1 用左除命令（m=A\b）来寻求矩阵m的最小二乘解
% coord_ini_matr=[];coord_x_att=[];
% for i=1:length(best_disratio)
% tmp=
% coord_ini_matr=[coord_ini_matr;]
% tmp=ones(length(best_disratio),1);
tmp=ones(size(best_disratio,1),1);
coord_ini_matr=[[Ipts1_impv(cor1_1).y]',[Ipts1_impv(cor1_1).x]',tmp];
coord_x_att=[Ipts2_impv(cor2_1).y]';
att_matr_13=coord_ini_matr\coord_x_att;
coord_y_att=[Ipts2_impv(cor2_1).x]';
att_matr_46=coord_ini_matr\coord_y_att;
att_matr_ls=[att_matr_13';att_matr_46';0,0,1];
% T_m{num,1}=att_matr_ls;
%%% 真实的放射攻击矩阵
%%%% 法2：先用已知X坐标拟合m11,m12,m13,再用已知Y坐标拟合m21,m22,m23,并和法1的结果相比较，验证矩阵左除是否是LS解
% att_matr_13_1=lsqnonneg(coord_ini_matr,coord_x_att);
% att_matr_46_1=lsqnonneg(coord_ini_matr,coord_y_att);
% att_matr_ls1=[att_matr_13_1';att_matr_46_1';0,0,1]; %%这不是最小二乘法吧？？？和att_matr_ls差别非常大
% Make vectors with the coordinates of the best matches
Pos1=[[Ipts1_impv(cor1_1).y]',[Ipts1_impv(cor1_1).x]'];
Pos2=[[Ipts2_impv(cor2_1).y]',[Ipts2_impv(cor2_1).x]'];
% Pos1=Pos1(1:15,:);
% Pos2=Pos2(1:15,:);
% Show both images
% I = zeros([size(I1,1) size(I1,2)*2 size(I1,3)]);
% I(:,1:size(I1,2),:)=I1; I(:,size(I1,2)+1:size(I1,2)+size(I2,2),:)=I2;
% figure, imshow(I); hold on;
%%gray image
I = zeros([round(1.6*size(I1,1)) round(size(I1,2)*2.6) ]);
I(1:size(I1,1),1:size(I1,2))=I1; I(1:size(I2,1),size(I1,2)+1:size(I1,2)+size(I2,2))=I2;
figure, imshow(I,[]); hold on;
% Show the best matches
plot([Pos1(:,2) Pos2(:,2)+size(I1,2)]',[Pos1(:,1) Pos2(:,1)]','-');
plot([Pos1(:,2) Pos2(:,2)+size(I1,2)]',[Pos1(:,1) Pos2(:,1)]','o');
%% 等比例缩放与旋转、不等比例缩放与旋转、仅有不等比例缩放的判定，若为不等比例缩放与旋转的攻击，则判定先R还是先S
m11=att_matr_ls(1,1);m12=att_matr_ls(1,2);
m21=att_matr_ls(2,1);m22=att_matr_ls(2,2);
% Sca_flag_un=0; %%不等比例缩放标志，对于先R在不等比例S，必须先还原S，然后再还原角度
% Sca_flag=0;
% RS_dist_condi1=abs(abs(m12/m11)-abs(m21/m22)); %% 到底是先不等比S还先R的区别条件
% % RS_dist_condi2=abs(abs(m11/m22)-abs(m21/m12)); %% 备注：此处在多想一下，有没有更精准的判定准则？？？
% % RS_dist_condi2 判定条件没多大意义，不等比例S、R的组合攻击都是如此
% if abs(m11-m22)>=0.05 %% 存在不等比例缩放或不等比例缩放与旋转的攻击
% if RS_dist_condi1<=0.015 %% 判定为先R再S（不等比）Plane这幅图比较特殊，这个值需大一些
% % [Sx,theta1]=solve('Sx*cosd(theta1)=m11','-Sx*sind(theta1)=m12','Sx','theta1');
% % [Sy,theta2]=solve('Sy*sind(theta2)=m21','Sy*sind(theta2)=m22','Sy','theta2');
% theta1=atand(-m12/m11); theta2=atand(m21/m22);
% Sx_nu=sqrt(m11^2+m12^2);Sy_nu=sqrt(m21^2+m22^2);
% Sx_nu
% Sy_nu
% theta=(theta1+theta2)/2;
% theta
% if abs(theta)<0.25
% disp('仅存在不等比例缩放');
% Sx_nu=m11;Sy_nu=m22;
% Sx_nu
% Sy_nu
% Rot_flag=0;
% Sca_flag=1; %%此处还是用Sca_flag吧，为了和水印判断那段程序统一
% [ro2,co2]=size(I2);
% x_s_rev=round(ro2/Sx_nu);y_s_rev=round(co2/Sy_nu);
% I4=imresize(im2uint8(I2),[x_s_rev y_s_rev],'bicubic'); %%将不等比例缩放还原
% Sx_rv1(num,1)=Sx_nu; Sy_rv1(num,1)=Sy_nu;
%
%
% else
% disp('存在不等比例缩放与旋转的攻击(先R再S)');
% [ro2,co2]=size(I2);
% x_s_rev=round(ro2/Sx_nu);y_s_rev=round(co2/Sy_nu); %%对于先R再S的不等比例攻击必须先还原S，若先还原角度则会造成类似于Shearing的攻击
% I2=imresize(im2uint8(I2),[x_s_rev y_s_rev],'bicubic'); %%备注：论文中如何描述、阐明？？？
% Sca_flag=1; %%此处还是先用Sca_flag吧，为了和水印判断那段程序统一
% I3_tmp=imrotate(im2uint8(I2),-theta,'bicubic','crop');
% I3=im2double(I3_tmp);
% figure,imshow(I3,[]),title('旋转theta角度后校正的图像');
% Rot_flag=1; %imrotate攻击标志
% I4=im2uint8(I3);
% Sx_rv1(num,1)=Sx_nu; Sy_rv1(num,1)=Sy_nu;
% thet_rv1(num,1)=theta;
% end
% else %% 判定为先S再R（不等比）
% % [Sx,theta1]=solve('Sx*cosd(theta1)=m11','Sx*sind(theta1)=m21','Sx','theta1');
% % [Sy,theta2]=solve('-Sy*sind(theta2)=m12','Sy*cosd(theta2)=m22','Sy','theta2');
% theta1=atand(m21/m11); theta2=atand(-m12/m22);
% Sx_nu=sqrt(m11^2+m21^2);Sy_nu=sqrt(m12^2+m22^2);
% Sx_nu
% Sy_nu
% theta=(theta1+theta2)/2;
% theta
% if abs(theta)<0.25
% disp('仅存在不等比例缩放');
% Sx_nu=m11;Sy_nu=m22;
% Sx_nu
% Sy_nu
% Rot_flag=0;
% Sca_flag=1; %%此处还是用Sca_flag吧，为了和水印判断那段程序统一
% [ro2,co2]=size(I2);
% x_s_rev=round(ro2/Sx_nu);y_s_rev=round(co2/Sy_nu);
% I4=imresize(im2uint8(I2),[x_s_rev y_s_rev],'bicubic'); %%将不等比例缩放还原
% Sx_rv1(num,1)=Sx_nu; Sy_rv1(num,1)=Sy_nu;
% else
% disp('存在不等比例缩放与旋转的攻击(先S再R)')
% I3_tmp=imrotate(im2uint8(I2),-theta,'bicubic','crop');
% I3=im2double(I3_tmp);
% figure,imshow(I3,[]),title('旋转theta角度后校正的图像'); %%%先还原角度
% Rot_flag=1; %imrotate攻击标志
% [ro3,co3]=size(I3);
% x_s_rev=round(ro3/Sx_nu);y_s_rev=round(co3/Sy_nu);
% I3=imresize(im2uint8(I3),[x_s_rev y_s_rev],'bicubic'); %%%再还原缩放尺度
% I4=I3;
% Sca_flag=1; %%此处还是先用Sca_flag吧，为了和水印判断那段程序统一
% Sx_rv1(num,1)=Sx_nu; Sy_rv1(num,1)=Sy_nu;
% thet_rv1(num,1)=theta;
% end
% end
% else %%等比例缩放或旋转或等比例缩放与旋转攻击同时有
%
%
% theta1=atand(-m12/m11); theta2=atand(m21/m22);theta3=atand(m21/m11);theta4=atand(-m12/m22);
% theta=(theta1+theta2+theta3+theta4)/4;
% theta
%
% S_un1=sqrt(m11^2+m12^2);S_un2=sqrt(m21^2+m22^2);S_un3=sqrt(m11^2+m21^2);S_un4=sqrt(m12^2+m22^2);
% S_un=(S_un1+S_un2+S_un3+S_un4)/4;
% S_un
%
% if abs(theta)<=0.25
% % figure,imshow(I2,[]);title('没受到旋转攻击的图片');
% I3=I2;
% Rot_flag=0;
% else
% I3_tmp=imrotate(im2uint8(I2),-theta,'bicubic','crop');
% I3=im2double(I3_tmp);
% figure,imshow(I3,[]),title('旋转theta角度后校正的图像');
% Rot_flag=1; %imrotate攻击标志
% thet_rv1(num,1)=theta;
% end
% %%判定有无等比例缩放攻击，若有则进行相应的还原
% if abs(S_un-1.0)<=0.015
% figure,imshow(I3);title('没受到等比例缩放攻击的图片');
% I4=im2uint8(I3);
% % Sca_flag=0;
% else
% [ro3,co3]=size(I3);
% x_s_rev=round(ro3/S_un);y_s_rev=round(co3/S_un);
% I4=imresize(im2uint8(I3),[x_s_rev y_s_rev],'bicubic');
% figure,imshow(uint8(I4),[ ]),title('缩放或旋转缩放攻击校正后的图像');
% Sca_flag=1; %受到缩放攻击的标志
% Sca_flag_un=1;
% S_un_rv1(num,1)=S_un;
% end
%
% end
%% 求出仿射矩阵T，利用求逆来还原图像
att_matr_ls1=[m22,m12,0;m21,m11,0;0,0,1];
T_m{num,1}=att_matr_ls1;
T_inv=inv(att_matr_ls1);
transa1=maketform('affine',T_inv);
I4=imtransform(im2uint8(I2),transa1);%WIa为遭遇RST几何攻击后的图像
figure;imshow(uint8(I4),[]);
[ro4,co4]=size(I4);
% I4=im2uint8(I4);
m=0;I4_1=[];
for i=1:ro4
if sum(uint8(I4(i,:)))>512*40
m=m+1;
I4_1(m,:)=I4(i,:);
end
end
n=0;
for i=1:co4
if sum(I4_1(:,i))>512*40
n=n+1;
I4_2(:,n)=I4_1(:,i);
end
end
I4=I4_2;
%% 首先在空域寻找嵌入特殊信息的位置块，以便定位到边缘被破坏的图片的边界以及再一次精准校正
Nc_total4=[];nc_coord_tal4=[];coord_avg4=[];
coordx=25;coordy=25;
[Nc_sta,nc_coord_alt,coord_x_y]=odd_even_ext24_66_1(I4,coordx,coordy,flg_grp(1,:));
Nc_total4=[Nc_total4;Nc_sta];nc_coord_tal4=[nc_coord_tal4;nc_coord_alt];coord_avg4=[coord_avg4;coord_x_y];
coordx=25;coordy=463;
[Nc_sta,nc_coord_alt,coord_x_y]=odd_even_ext24_lena12(I4,coordx,coordy,flg_grp(2,:));
Nc_total4=[Nc_total4;Nc_sta];nc_coord_tal4=[nc_coord_tal4;nc_coord_alt];coord_avg4=[coord_avg4;coord_x_y];
% coordx=247;coordy=247;
% [Nc_sta,nc_coord_alt,coord_x_y]=odd_even_ext24_66(I5,coordx,coordy,flg_grp(3,:));
% Nc_total=[Nc_total;Nc_sta];nc_coord_tal=[nc_coord_tal;nc_coord_alt];coord_avg=[coord_avg;coord_x_y];
coordx=463;coordy=25;
[Nc_sta,nc_coord_alt,coord_x_y]=odd_even_ext24_lena21(I4,coordx,coordy,flg_grp(4,:));
Nc_total4=[Nc_total4;Nc_sta];nc_coord_tal4=[nc_coord_tal4;nc_coord_alt];coord_avg4=[coord_avg4;coord_x_y];
coordx=463;coordy=463;
[Nc_sta,nc_coord_alt,coord_x_y]=odd_even_ext24_66_1(I4,coordx,coordy,flg_grp(5,:));
Nc_total4=[Nc_total4;Nc_sta];nc_coord_tal4=[nc_coord_tal4;nc_coord_alt];coord_avg4=[coord_avg4;coord_x_y];
%%利用特殊位置的点相对距离求出像素的偏差，并且再一次精准校正
x1x4_dis=coord_avg4(3,1)-coord_avg4(1,1);
x2x5_dis=coord_avg4(4,1)-coord_avg4(2,1);
x_devi=(438-(x1x4_dis+x2x5_dis)*0.5)/(438/512); %若x_devi为正则需要放大，为负则需要缩小
y1y4_dis=coord_avg4(2,2)-coord_avg4(1,2);
y2y5_dis=coord_avg4(4,2)-coord_avg4(3,2);
y_devi=(438-(y1y4_dis+y2y5_dis)*0.5)/(438/512); %若y_devi为正则需要放大，为负则需要缩小
[ro4,co4]=size(I4);
if abs(x_devi)>=2.0&&abs(y_devi)>=2.0
ro5=round(x_devi)+ro4;co5=round(y_devi)+co4;
I5=imresize(I4,[ro5 co5],'bicubic');
Scl_fc(count,1)=1;
elseif abs(x_devi)>=2.0&&abs(y_devi)<2.0
3 运行结果
4 参考文献
[1]尤鸿霞, 郭江峰. 基于Matlab的二值图象数据隐藏算法的实现[J]. 南通纺织职业技术学院学报, 2008.
[2]庞镇明. 基于MATLAB小波变换在图像处理中算法的实现[D]. 北京理工大学.
