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 運作結果
