% 由高斯正態分布序列 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('攻擊后的水印圖像');