图像数字水印Matlab源代码

1、运行MAIN.m即可开始水印的嵌入和提取。
2、文件夹中的两幅图片为载体图像和水印图像。
3、其他文件为主程序所调用的自定义函数，说明如下：
sdwt.m:对图像依视觉能量进行树状小波分解
embed.m:对标记的嵌入点进行水印嵌入
nembed.m:对每个节点实施水印嵌入
sidwt.m:对嵌入后的树形子图以小波逆变换进行重组
sdwt_ex.m:依密钥树对含水印图像进行分解
extract.m:依密钥树抽取水印
nextract.m:对每个节点实施水印抽取
jadeR.m:JADE算法，用于实现ICA
fuse_pca.m:PCA算法，用于实现融合
rand_orth.m:生成混叠矩阵随机数





MAIN.m 主程序

%-------------------水印嵌入------------------------------------------------
while 1
clear;
c=0.3;
a=imread('lina.jpg');%原图像
b=imread('changsha.bmp')*255;%二值水印图像
[m1,n1]=size(a);
[m2,n2]=size(b);
e0=(sum(sum(a.^2)))/(m1*n1);
e0=c*e0;%计算基准能量
[t,tkey]=sdwt(double(a),'db2',m2,n2,e0);%树状小波分解
[t,tkey]=embed(t,tkey,b);%嵌入水印
aw=sidwt(t,'db2');%重组
figure(1);
subplot(1,2,1);imshow(uint8(a));title('原图');
subplot(1,2,2);imshow(uint8(aw));title('嵌入图');
imwrite(uint8(aw),'watermark.jpg');
% csvwrite('key.txt',reshape(tkey,m2,n2));
v1=m1*m1*255*255;
v2=sum(sum((double(a)-aw).^2));
snr=10*log10(v1/v2);% 峰值信噪比snr。
disp('信噪比');disp(snr);
%---攻击预处理--------------------------------------------------------------
% clear;
% pause;
f=imread('watermark.jpg');
%将含水印图像f归一化,以便于攻击处理。
m=max(max(f));
f=double(f)./double(m);
%---攻击-------------------------------------------------------------------
attack=0;
switch attack
case 0,
attackf=f;
att='未攻击';
case 1,
%%1. JPEG 压缩
imwrite(f,'attackf.jpg','jpg','quality',30);
attackf=imread('attackf.jpg');
attackf=double(attackf)/255;
att='JPEG压缩';
case 2,
% %2. 高斯低通滤波
h=fspecial('gaussian',3,1);
attackf=filter2(h,f);
att='高斯低通滤波';
case 3,
%%3. 直方图均衡化
attackf=histeq(f);
att='直方图均衡化';
case 4,
%%4. 图像增亮
attackf=imadjust(f,[],[0.4,1]);
att='图像增亮';
case 5,
%%5. 图像变暗
attackf=imadjust(f,[],[0,0.85]);
att='图像变暗';
case 6,
%%6. 增加对比度
attackf=imadjust(f,[0.3,0.6],[]);
att='增加对比度';
case 7,
%%7. 降低对比度
attackf=imadjust(f,[],[0.2,0.8]);
att='降低对比度';
case 8,
%%8. 添加高斯噪声
attackf=imnoise(f,'gaussian',0,0.01);
att='添加高斯噪声';
case 9,
%%9. 椒盐噪声
attackf=imnoise(f,'salt & pepper',0.06);
att='椒盐噪声';
case 10,
%%10. 添加乘积性噪声
attackf=imnoise(f,'speckle',0.08);
att='添加乘积性噪声';
end;
%---攻击后处理--------------------------------------------------------------
f=attackf.*double(m);
figure(2);
imshow(uint8(f));%显示水印嵌入图攻击后效果
title(att);
imwrite(uint8(f),'watermark.jpg');
%---提取水印----------------------------------------------------------------
% clear;
a=imread('watermark.jpg');
t=sdwt_ex(double(a),'db2',tkey);%根据密钥树分解
[w,map]=extract(t,tkey);%抽取水印
[r,c]=size(w);
figure(3);
for i=1:r
subplot(ceil(r/3),3,i)
imshow(255-100*abs(uint8(reshape(w(i,:),map(1),map(2)))));
title(strcat('抽取水印图',num2str(i)));
end;
%---分析水印并融合----------------------------------------------------------
Rarr=[];
for i=1:r %建立相关系数矩阵
for j=i+1:r
R(i,j)=abs(corr2(reshape(w(i,:),map(1),map(2)),reshape(w(j,:),map(1),map(2))));
if i~=j
Rarr=[Rarr;i,j,abs(R(i,j))];
end;
end;
end;

[cor,idx]=sort(Rarr(:,3),'descend');
max1=cor(1);max2=cor(2);
m1=Rarr(idx(1),1);n1=Rarr(idx(1),2);
m2=Rarr(idx(2),1);n2=Rarr(idx(2),2);
maxcor(1)=abs(corr2(reshape(w(m1,:),map(1),map(2)),reshape(w(m2,:),map(1),map(2))))
maxcor(2)=abs(corr2(reshape(w(m1,:),map(1),map(2)),reshape(w(n2,:),map(1),map(2))))
maxcor(3)=abs(corr2(reshape(w(n1,:),map(1),map(2)),reshape(w(m2,:),map(1),map(2))))
maxcor(4)=abs(corr2(reshape(w(n1,:),map(1),map(2)),reshape(w(n2,:),map(1),map(2))))
if mean(maxcor)>(max1*0.8)
stdw=fuse_pca(reshape(w(m1,:),map(1),map(2)), reshape(w(n1,:),map(1),map(2)));
figure(4);
subplot(1,2,1)
imshow(b);
title('原水印');
subplot(1,2,2)
wf=255-100*abs(uint8(stdw));
imshow(wf);
title('pca融合后的水印');
else
disp('重新融合');
end;
if abs(corr2(wf,b))>0.9
break;
end;
end;
disp('水印相似度');disp(corr2(wf,b));
%--------------------------------------------------------------------------

sdwt.m:对图像依视觉能量进行树状小波分解

function [t,tkey]=sdwt(node,wname,m2,n2,e0)
[A,H,V,D]=dwt2(node,wname);
[mi,ni]=size(A);
if mi*ni<m2*n2
t=node;tkey=1;
return;
end;
ea=(sum(sum(A.^2)))/(mi*ni);
eh=(sum(sum(H.^2)))/(mi*ni);
ev=(sum(sum(V.^2)))/(mi*ni);
ed=(sum(sum(D.^2)))/(mi*ni);
if ea<e0&ea>0.1*e0
[A,ka]=sdwt(A,wname,m2,n2,e0);
else
ka=-1;
end;
if eh<e0&eh>0.1*e0
[H,kh]=sdwt(H,wname,m2,n2,e0);
else
kh=-1;
end;
if ev>e0&ev>0.1*e0
[V,kv]=sdwt(V,wname,m2,n2,e0);
else
kv=-1;
end;
if ed<e0&ed>0.1*e0
[D,kd]=sdwt(D,wname,m2,n2,e0);
else
kd=-1;
end;
t={A,H,V,D};
tkey={ka,kh,kv,kd};

return;

embed.m:对标记的嵌入点进行水印嵌入

function [t,tkey]=embed(t,tkey,b)
[m2,n2]=size(b);
for i=1:4
if ~iscell(tkey{1,i})
if tkey{1,i}==1
[t{1,i},tkey{1,i}]=nembed(t{1,i},b);
end;
else
[t{1,i},tkey{1,i}]=embed(t{1,i},tkey{1,i},b);
end;
end;
return;

nembed.m:对每个节点实施水印嵌入

function [wnode,nodekey]=nembed(nd,b)
[m1,n1]=size(nd);
[m2,n2]=size(b);
node=nd;
node1=node(1:m2*n2);
node2=node(m2*n2+1:m1*n1);
b=b(1:m2*n2);

%---互信息法---------------------------------------------------------------
S=[];
for i=1:2
switch i
case 1, s=double(node1);
case 2, s=double(b);
end
s=s-mean(s); % 归一
s=s/std(s,1); % 标准化
S=[S; s];
end
A=rand_orth(2);
% A=rand(2,2);
X=A*S;
x1=X(1,:);
x1=[x1,node2];
x2=X(2,:)
nodekey=reshape(x2,m2,n2);
%--------------------------------------------------------------------------
%---加法-------------------------------------------------------------------
% x1=node1+b*0.02;
% x1=[x1 node2];
% nodekey=[0.02,-1];
%--------------------------------------------------------------------------
wnode=reshape(x1,m1,n1);

return;

sidwt.m:对嵌入后的树形子图以小波逆变换进行重组
function t=sidwt(t,wname)
if ~iscell(t{1,1})
if ~iscell(t{1,2})
if ~iscell(t{1,3})
if ~iscell(t{1,4})
[mi(1),ni(1)]=size(t{1,1});[mi(2),ni(2)]=size(t{1,2});
[mi(3),ni(3)]=size(t{1,3});[mi(4),ni(4)]=size(t{1,4});
m=min(mi);n=min(ni);
t{1,1}=t{1,1}(1:m,1:n);
t{1,2}=t{1,2}(1:m,1:n);
t{1,3}=t{1,3}(1:m,1:n);
t{1,4}=t{1,4}(1:m,1:n);
t=idwt2(t{1,1},t{1,2},t{1,3},t{1,4},wname);
else
t{1,4}=sidwt(t{1,4},wname);
[mi(1),ni(1)]=size(t{1,1});[mi(2),ni(2)]=size(t{1,2});
[mi(3),ni(3)]=size(t{1,3});[mi(4),ni(4)]=size(t{1,4});
m=min(mi);n=min(ni);
t{1,1}=t{1,1}(1:m,1:n);
t{1,2}=t{1,2}(1:m,1:n);
t{1,3}=t{1,3}(1:m,1:n);
t{1,4}=t{1,4}(1:m,1:n);
t=idwt2(t{1,1},t{1,2},t{1,3},t{1,4},wname);
end;
else
if iscell(t{1,4})
t{1,4}=sidwt(t{1,4},wname);
end;
t{1,3}=sidwt(t{1,3},wname);
[mi(1),ni(1)]=size(t{1,1});[mi(2),ni(2)]=size(t{1,2});
[mi(3),ni(3)]=size(t{1,3});[mi(4),ni(4)]=size(t{1,4});
m=min(mi);n=min(ni);
t{1,1}=t{1,1}(1:m,1:n);
t{1,2}=t{1,2}(1:m,1:n);
t{1,3}=t{1,3}(1:m,1:n);
t{1,4}=t{1,4}(1:m,1:n);
t=idwt2(t{1,1},t{1,2},t{1,3},t{1,4},wname);

end;
else
if iscell(t{1,4})
t{1,4}=sidwt(t{1,4},wname);
end;
if iscell(t{1,3})
t{1,3}=sidwt(t{1,3},wname);
end;
t{1,2}=sidwt(t{1,2},wname);
[mi(1),ni(1)]=size(t{1,1});[mi(2),ni(2)]=size(t{1,2});
[mi(3),ni(3)]=size(t{1,3});[mi(4),ni(4)]=size(t{1,4});
m=min(mi);n=min(ni);
t{1,1}=t{1,1}(1:m,1:n);
t{1,2}=t{1,2}(1:m,1:n);
t{1,3}=t{1,3}(1:m,1:n);
t{1,4}=t{1,4}(1:m,1:n);
t=idwt2(t{1,1},t{1,2},t{1,3},t{1,4},wname);
end;
else
if iscell(t{1,4})
t{1,4}=sidwt(t{1,4},wname);
end;
if iscell(t{1,3})
t{1,3}=sidwt(t{1,3},wname);
end;
if iscell(t{1,2})
t{1,2}=sidwt(t{1,2},wname);
end;
t{1,1}=sidwt(t{1,1},wname);
[mi(1),ni(1)]=size(t{1,1});[mi(2),ni(2)]=size(t{1,2});
[mi(3),ni(3)]=size(t{1,3});[mi(4),ni(4)]=size(t{1,4});
m=min(mi);n=min(ni);
t{1,1}=t{1,1}(1:m,1:n);
t{1,2}=t{1,2}(1:m,1:n);
t{1,3}=t{1,3}(1:m,1:n);
t{1,4}=t{1,4}(1:m,1:n);
t=idwt2(t{1,1},t{1,2},t{1,3},t{1,4},wname);
end;
return;

sdwt_ex.m:依密钥树对含水印图像进行分解
function t=sdwt_ex(node,wname,tkey)
[A,H,V,D]=dwt2(node,wname);
AHVD={A,H,V,D};
for i=1:length(tkey)
if iscell(tkey{1,i})
t{1,i}=sdwt_ex(AHVD{1,i},wname,tkey{1,i});
else
if tkey{1,i}~=-1
t{1,i}=node;
else
t{1,i}={};
end;
end;
end;

extract.m:依密钥树抽取水印

function [w,map]=extract(t,tkey);
w=[];
for i=1:4
if ~iscell(tkey{1,i})
if tkey{1,i}~=-1
[wi,map]=nextract(t{1,i},tkey{1,i});
else
wi=[];
end;
else
[wi,map]=extract(t{1,i},tkey{1,i});
end;
if length(wi)>0
w=[w;wi];
end;
end;
return;

nextract.m:对每个节点实施水印抽取

function [wi,map]=nextract(tnd,tkey)
[m2,n2]=size(tkey);
x1=tnd(1:m2*n2);
x2=tkey(1:m2*n2);
X=[x1;x2];
B=jadeR(double(X),2);
Y=B*double(X);
y2=Y(2,:);
wi=y2;
map=[m2,n2];
return;

fuse_pca.m:PCA算法，用于实现融合

function Y = fuse_pca(M1, M2)
[z1 s1] = size(M1);
[z2 s2] = size(M2);
if (z1 ~= z2) | (s1 ~= s2)
error('Input images are not of same size');
end;
[V, D] = eig(cov([M1(:) M2(:)]));
if (D(1,1) > D(2,2))
a = V(:,1)./sum(V(:,1));
else
a = V(:,2)./sum(V(:,2));
end;
Y = a(1)*M1+a(2)*M2;

rand_orth.m:生成混叠矩阵随机数

function W=rand_orth(n,m);
if (nargin<2)
m=n;
end
W=rand(m)-.5;
[W,cococococo]=qr(W);
W=W(1:m,1:n);

jadeR.m:JADE算法，用于实现ICA

function B = jadeR(X,m)
% B = jadeR(X, m) is an m*n matrix such that Y=B*X are separated sources
% extracted from the n*T data matrix X.
% If m is omitted, B=jadeR(X) is a square n*n matrix (as many sources as sensors)
%
% Blind separation of real signals with JADE. Version 1.8. * If X is an nxT data matrix (n sensors, T samples) then
% B=jadeR(X) is a nxn separating matrix such that S=B*X is an nxT
% matrix of estimated source signals.
% * If B=jadeR(X,m), then B has size mxn so that only m sources are
% extracted. This is done by restricting the operation of jadeR
% to the m first principal components.
% * Also, the rows of B are ordered such that the columns of pinv(B)
% are in order of decreasing norm; this has the effect that the
% `most energetically significant' components appear first in the
% rows of S=B*X.
%
% Quick notes (more at the end of this file)
%
% o this code is for REAL-valued signals. An implementation of JADE
% for both real and complex signals is also available from
% http://sig.enst.fr/~cardoso/stuff.html
%
% o This algorithm differs from the first released implementations of
% JADE in that it has been optimized to deal more efficiently
% 1) with real signals (as opposed to complex)
% 2) with the case when the ICA model does not necessarily hold.
%
% o There is a practical limit to the number of independent
% components that can be extracted with this implementation. Note
% that the first step of JADE amounts to a PCA with dimensionality
% reduction from n to m (which defaults to n). In practice m
% cannot be `very large' (more than 40, 50, 60... depending on
% available memory)
%
% o See more notes, references and revision history at the end of
% this file and more stuff on the WEB
% http://sig.enst.fr/~cardoso/stuff.html
%
% o This code is supposed to do a good job! Please report any
% problem to cardoso@sig.enst.fr

cardoso@sig.enst.fr

verbose = 1 ; % Set to 0 for quiet operation

% Finding the number of sources
[n,T] = size(X);
if nargin==1, m=n ; end; % Number of sources defaults to # of sensors
if m>n , fprintf('jade -> Do not ask more sources than sensors here!!!\n'), return,end
if verbose, fprintf('jade -> Looking for %d sources\n',m); end ;

% to do: add a warning about complex signals

% Mean removal
%=============
if verbose, fprintf('jade -> Removing the mean value\n'); end
X = X - mean(X')' * ones(1,T);

%%% whitening & projection onto signal subspace
% ===========================================
if verbose, fprintf('jade -> Whitening the data\n'); end

[U,D] = eig((X*X')/T) ; %% An eigen basis for the sample covariance matrix
[Ds,k] = sort(diag(D)) ; %% Sort by increasing variances
PCs = ; %% The m most significant princip. comp. by decreasing variance

%% --- PCA ----------------------------------------------------------
B = ; % At this stage, B does the PCA on m components

%% --- Scaling ------------------------------------------------------
scales = sqrt(Ds(PCs)) ; % The scales of the principal components .
B = diag(1./scales)*B ; % Now, B does PCA followed by a rescaling = sphering

%% --- Sphering ------------------------------------------------------
X = B*X;

Any further rotation of X will preserve the
%%% property that X is a vector of uncorrelated components. It remains to find the
%%% rotation matrix such that the entries of X are not only uncorrelated but also `as
%%% independent as possible'. This independence is measured by correlations of order
%%% higher than 2. We have defined such a measure of independence which
%%% 1) is a reasonable approximation of the mutual information
%%% 2) can be optimized by a `fast algorithm'
%%% This measure of independence also corresponds to the `diagonality' of a set of
%%% cumulant matrices. The code below finds the `missing rotation ' as the matrix which
%%% best diagonalizes a particular set of cumulant matrices.

%%% Estimation of the cumulant matrices.
% ====================================
if verbose, fprintf('jade -> Estimating cumulant matrices\n'); end

%% Reshaping of the data, hoping to speed up things a little bit...
X = X';

dimsymm = (m*(m+1))/2; % Dim. of the space of real symm matrices
nbcm = dimsymm ; % number of cumulant matrices
CM = zeros(m,m*nbcm); % Storage for cumulant matrices
R = eye(m); %%
Qij = zeros(m); % Temp for a cum. matrix
Xim = zeros(m,1); % Temp
Xijm = zeros(m,1); % Temp
Uns = ones(1,m); % for convenience

%% I am using a symmetry trick to save storage. I should write a short note one of these
%% days explaining what is going on here.
%%
Range = % will index the columns of CM where to store the cum. mats.

for im = Xim =
Xijm= Xim.*Xim ;
%% the joint diagonalization criterion
Qij =
= Qij ;
Range = Range + m ;
for jm = Xijm =
Qij =
= Qij ;
Range = Range + m ;
end ;
end;
%%%% Now we have nbcm = It can be found at
%% "http://www.tsi.enst.fr/~cardoso/Papers.PS/neuralcomp_2ppf.ps",
%%
%%
%%
%% if 1,
%%
%% Matcum = zeros(m,m,m,m) ;
%% for i1=1:m,
%% for i2=1:m,
%% for i3=1:m,
%% for i4=1:m,
%% Matcum(i1,i2,i3,i4) = - R(i1,i2)*R(i3,i4) ...
%% - R(i1,i3)*R(i2,i4) ...
%% - R(i1,i4)*R(i2,i3) ;
%% end
%% end
%% end
%% end
%%
%% %% Step 2; We compute a basis of the space of symmetric m*m matrices
%% CMM = zeros(m, m, nbcm) ; %% Holds the basis.
%% icm = 0 ; %% index to the elements of the basis
%% vi = zeros(m,1); %% the ith basis vetor of R^m
%% vj = zeros(m,1); %% the jth basis vetor of R^m
%% Id = eye (m) ; %% convenience
%% for im=1:m,
%% vi =
%% icm = icm + 1 ;
%% CMM(:, :, icm) = vi*vi' ;
%% for jm=1:im-1,
%% vj =
%% icm = icm + 1 ;
%% CMM(:, :, icm) = sqrt(0.5) * (vi*vj'+vj*vi') ;
%% end
%% end
%%
%% %% Step 3. We compute the image of each basis element by the cumulant tensor and store it back into CMM.
%% mat = zeros(m) ; %% tmp
%% for icm=1:nbcm
%% mat =
%% for i1=1:m
%% for i2=1:m
%% CMM(i1, i2, icm) = .* mat )) ;
%% end
%% end
%% end;
%% %% This is doing something like \sum_kl [ Cum(xi,xj,xk,xl) * mat_kl ]
%%
%% %% Step 4. Now, we can check that CMM and CM are equivalent
%% Range =
%% for icm=1:nbcm,
%% M1 =
%% M2 =
%% Range = Range + m ;
%% norm (M1-M2, 'fro' ) , %% This should be a numerical zero.
%% end;
%%
%% end; %% End of the demo code for the computation of cumulant matrices
%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%