水印源码.txt

用matlab编写的产生水印以及若干攻击测试的代码

% 由高斯正态分布序列 g1 产生 36×4 的水印信
%号 w0，w0 由(0,1)组成。
clear
randn('state',1106);
g1=randn(36,4);
for i=1:36
     for j=1:4
          if g1(i,j)>=0
                w0(i,j)=1;
         else w0(i,j)=0;
         end;
     end;
end;
figure;
imshow(w0);title('水印');


% 对水印信号 w0 进行(7,4)汉明编码,得到一
%36×7 的分组码 x0。
x0=w0;
for i=1:36
    s=8*x0(i,1)+4*x0(i,2)+2*x0(i,3)+x0(i,4);
     switch s
         case 0
              x0(i,5)=0;x0(i,6)=0;x0(i,7)=0;
         case 1
              x0(i,5)=0;x0(i,6)=1;x0(i,7)=1;
         case 2
              x0(i,5)=1;x0(i,6)=1;x0(i,7)=0;
         case 3
              x0(i,5)=1;x0(i,6)=0;x0(i,7)=1;
         case 4
              x0(i,5)=1;x0(i,6)=1;x0(i,7)=1;
         case 5
              x0(i,5)=1;x0(i,6)=0;x0(i,7)=0;
         case 6
              x0(i,5)=0;x0(i,6)=0;x0(i,7)=1;
         case 7
              x0(i,5)=0;x0(i,6)=1;x0(i,7)=0;
         case 8
              x0(i,5)=1;x0(i,6)=0;x0(i,7)=1;
         case 9
              x0(i,5)=1;x0(i,6)=1;x0(i,7)=0;
         case 10
              x0(i,5)=0;x0(i,6)=1;x0(i,7)=1;
         case 11
              x0(i,5)=0;x0(i,6)=0;x0(i,7)=0;
         case 12
              x0(i,5)=0;x0(i,6)=1;x0(i,7)=0;
         case 13
              x0(i,5)=0;x0(i,6)=0;x0(i,7)=1;
         case 14
              x0(i,5)=1;x0(i,6)=0;x0(i,7)=0;
         case 15
              x0(i,5)=1;x0(i,6)=1;x0(i,7)=1;
     end;
end;
% 对 x0 进行行向位扩展，得到一个由(-1,1)组成
%的扩展序列 y。cr 为扩展因子。
cr=256;
for i=1:252
     if x0(i)==1
         y(i,1:cr)=1;
     else y(i,1:cr)=-1;
     end;
end;
y(253:256,:)=0;
% 以下产生伪随机序列 p。为此先设定密钥(1114)
%并产生高斯正态分布序列g2，再由g2产生由(-1,1)
%组成的伪随机序列 p。
randn('state',1114);
g2=randn(256,256);
for i=1:256
     for j=1:256
          if g2(i,j)>0
               p(i,j)=1;
          else p(i,j)=-1;
          end;
     end;
end;
yp=y.*p;
% 下面设定的每类块基准噪声阈值 jnd1 是通过
%反复实验确定出来的。
t1=1.7; t2=2.1; t3=2.5; t4=2.9;
t5=3.3; t6=3.7; t7=4.1; t8=4.5;
% 读入原图象并转换成双精度，k 是 8×8 图像
%块数。
f0=imread('cameraman.tif');
f0=double(f0);
[c,s]=size(f0);
k=c*s/64;
% 计算每块的方差 std2 和能量 e，并对能量 e 按
%5升序排列。
cf0=im2col(f0,[8,8],'distinct');
std1=std(cf0);
std2=std1.^2;
for i=1:k
e(i)=sum(cf0(:,i).^2);
end;
[e1,ind1]=sort(e);
% 按能量 e 的索引 ind1 顺序重排方差 std2，将结
%果存于 std3。
j=1;
for i=1:k
     z=ind1(i);
     std3(j)=std2(z);
     j=j+1;
end;
% 设定分类界限。
m1=median(std3(1:k/4));
m2=median(std3(k/4+1:k/2));
m3=median(std3(k/2+1:3*k/4));
m4=median(std3(3*k/4+1:k));
n1=e1(k/4);          n2=e1(2*k/4);
n3=e1(3*k/4);        n4=e1(k);
% 按能量 e 和方差 std2 将原图像块分成 8 类，并
%给每类块赋基准噪声阈值。
for i=1:k
     if e(i)<=n1
         if std2(i)>=m1
               jnd1(i)=t1;
         else jnd1(i)=t2;
         end;
    elseif (n1<e(i)&e(i)<=n2)
           if std2(i)>=m2
                jnd1(i)=t3;
           else jnd1(i)=t4;
           end;
     elseif(n2<e(i)&e(i)<=n3)
           if std2(i)>=m3
                jnd1(i)=t5;
           else jnd1(i)=t6;
           end;
     elseif (n3<e(i)&e(i)<=n4)
           if std2(i)>=m4
                jnd1(i)=t7;
           else jnd1(i)=t8;
           end;
     end;
end;
% 计算每类块的噪声阈值 jnd，它等于基准噪声
%阈值 jnd1 和附加噪声阈值 jnd2 之和。
deta=0.0035;
jnd2=deta*cf0;
for i=1:k
     jnd(:,i)=jnd1(i)+jnd2(:,i);
end;

% 嵌入已调制的水印信号 yp。并重构成含水印
%的图像 f1。
recf0=reshape(cf0,256,256);
rejnd=reshape(jnd,256,256);
recf1=recf0+rejnd.*yp;
% 重构嵌入水印的图像 f1。
cf1=reshape(recf1,64,1024);
f1=col2im(cf1,[8,8],[256,256],'distinct');
% 从含水印的图像 f1 中提取出 36×7 的分组码
%x1。
for i=1:252
     sk(i)=sum((recf1(i,:)-recf0(i,:)).*p(i,:));
end;
for i=1:252
     if sign(sk(i))==-1
           rex1(i)=0;
     else rex1(i)=1;
      end;
end;
x1=reshape(rex1,36,7);
% 对提取出来的分组码 x1 进行纠错解码，最后
%得到一个 36×4 的水印 w1。
for i=1:36
      s1(i)=x1(i,1)+x1(i,2)+x1(i,3)+x1(i,5);
      s2(i)=x1(i,2)+x1(i,3)+x1(i,4)+x1(i,6);
      s3(i)=x1(i,1)+x1(i,2)+x1(i,4)+x1(i,7);
      s1(i)=mod(s1(i),2);
      s2(i)=mod(s2(i),2);
      s3(i)=mod(s3(i),2);
      if (s1(i)==0&s2(i)==0)&(s3(i)==1)
           x1(i,7)=~x1(i,7);
      elseif (s1(i)==0&s2(i)==1)&(s3(i)==0)
           x1(i,6)=~x1(i,6);
      elseif (s1(i)==0&s2(i)==1)&(s3(i)==1)
           x1(i,4)=~x1(i,4);
      elseif (s1(i)==1&s2(i)==0)&(s3(i)==0)
           x1(i,5)=~x1(i,5);
      elseif (s1(i)==1&s2(i)==0)&(s3(i)==1)
           x1(i,1)=~x1(i,1);
      elseif (s1(i)==0&s2(i)==1)&(s3(i)==0)
           x1(i,3)=~x1(i,3);
      elseif (s1(i)==1&s2(i)==1)&(s3(i)==1)
           x1(i,2)=~x1(i,2);
      end;
end;
w1=x1(:,1:4);
% 计算原图像 f0 和嵌入水印的图像 f1 的信噪比
%snr。
v1=sum(sum(f0.^2));
v2=sum(sum((f0-f1).^2));
snr=10*log10(v1/v2);
% 计算原水印 w0 和提取出来的水印 w1 的相关
%系数 r。
r=corr2(w0,w1);
% 显示原始图像 f0 和嵌入水印的图像 f1，结果
%示于图 1、图 2。
figure
imshow(f0,[]);
title('原图像');
figure
imshow(f1,[]);
title('嵌入水印的图像');
% 将含水印图像 f1 归一化,以便于攻击处理。
m=max(max(f1));
f=f1/m;
%1. JPEG 压缩。
imwrite(f,'attackf.jpg','jpg','quality',30);
attackf=imread('attackf.jpg');
attackf=double(attackf)/255;
%2. 高斯低通滤波。
h=fspecial('gaussian',3,1);
attackf=filter2(h,f);
%3. 直方图均衡化。
attackf=histeq(f);
%4. 图像增亮。
attackf=imadjust(f,[],[0.4,1]);
%5. 图像变暗。
attackf=imadjust(f,[],[0,0.85]);
%6. 增加对比度。
attackf=imadjust(f,[0.3,0.6],[]);
%7. 降低对比度。
attackf=imadjust(f,[],[0.2,0.8]);
%8. 添加高斯噪声。
attackf=imnoise(f,'gaussian',0,0.01);
%9. 增加黑白像素点。
attackf=imnoise(f,'salt & pepper',0.06);
%10 添加乘积性噪声。
attackf=imnoise(f,'speckle',0.08);
% 显示攻击后的水印图像 attackf，为了简单起见
%这里只给出最后一种攻击即添加乘性噪声后的
%水印图像，其它攻击情况相同。
figure
imshow(attackf,[]);
title('攻击后的水印图像');