%基于HVS的混沌陣列水印算法
clear all
[I,map]=imread('d:\matpro\tuku\fds.bmp');
I=double(I);
%MuI=mean2(I)
%STDI=std2(I)
%En1=entropy(I)
%I=rgb2gray(I);
%map=rgb2gray(map);
%I=imread('rice.tif');
%load woman
invariancematrix1= abs(invariance(I,1,0));   %原始圖像的不變矩
I=double(I)/255;%歸一化
alpha=0.2;%嵌入水印的強度系數(shù)
dctm=dctmtx(8);
Fk=blkproc(I,[8 8],'P1*x*P2',dctm,dctm');
mask=[1 1 1 1 0 0 1 1 
      1 1 1 0 0 0 1 1 
      1 1 0 0 0 1 1 1
      1 0 0 0 1 1 1 1
      0 0 0 1 1 1 1 1
      0 0 1 1 1 1 1 1
      1 1 1 1 1 1 1 1
      1 1 1 1 1 1 1 1];
Fkc=blkproc(Fk,[8 8],'P1.*x',mask);
Cz=blkproc(Fk,[8 8],'chaoswm(x)');%提取次低頻DCT128x128系數(shù)陣 
w=choasq(0.200000000000101,128);%混沌水印的生成10E-14
%w=msqu(127,7);%M序列水印的生成
%Muw=mean2(w)
%STDw=std2(w)
%--------------------------------------
Cd=Cz+alpha.*abs(Cz).*w;%水印算法的乘法準則，v' = v + a x w; 0.01<a<100.
%Cd=Cz+alpha.*w;%加法準則
%--------------------------------------
Fkz=blkproc(Cd,[4 4],'chaoswn(x)'); 
Fki=Fkc+Fkz;
%------------------------------------------------
JND=hvsc(Fki);%視覺掩蓋的JND矩陣
Fkm=Fk.*(1-JND)+Fki.*JND;%視覺掩蓋的準則
%------------------------------------------------
I2=blkproc(Fkm,[8 8],'P1*x*P2',dctm',dctm);%圖像反變換
%I2=blkproc(Fki,[8 8],'P1*x*P2',dctm',dctm);
I2=(I2-min(min(I2)))./(max(max(I2))-min(min(I2)));%歸一化
invariancematrix2= abs(invariance(I2*255,1,0));   %原始圖像的不變矩
%------------------------------------------------
%HVS masking
%------------------------------------------------
%hvsa=hvsc(I2);%varying coffience
%hvsa=0.2;%constant coffience
%Im=I.*(1-hvsa)+I2.*hvsa;
%------------------------------------------------
%MSE=sqrt(sum(sum((I2*255-I*255).^2))/prod(size(I)));%均方根誤差(I'-I)^2/N^2;
%PSNR=20*log10(255/MSE);%水印算法中的峰峰信噪比，PSNR=10log(255^2/MSE^2),35.5125
%PSNR
%Ec=corr2(I,I2)
%------------------------------------------------
%imwrite(I2,'d:\matpro\chaos\c2.jpg','Quality',5);
imwrite(I2*255,map,'d:\matpro\chaos\c1.bmp');%圖像的存儲
%imwrite(I2,'d:\stf\matpro\chaos\c1.tif','tif');
%En2=entropy(I2)
%figure,imshow(I),figure,imshow(I2);
%------------------------------------------------
%水印的檢測過程
%------------------------------------------------
I3=imread('d:\matpro\chaos\c1.bmp');
%I3=double(I3);
Mud=mean2(I3);
STDd=std2(I3);
%I3=imread('cameraman.tif');
%I3=double(I3)/255;
%-------------------------------------------------
%figure,imshow(I3);
%en2=entropy(I3);
%-------------------------------------------------------------
%圖像的各種處理：加噪聲、低通濾波、中值濾波、剪切、比例縮放...
%-------------------------------------------------------------
%I3=imnoise(I3,'gaussian',0,0.07);
%I3=imnoise(I3,'salt&pepper',0.15);
%I3=imnoise(I3,'speckle',0.39);
%I3=wiener2(I3,[5 5]);
%I3=medfilt2(I3,[3 3]);
%J=dither(I3);
%I3=double(J);
%figure,imshow(J)
%------------------------------------------------
%image shearing
%------------------------------------------------
%I3=[I3(1:128,1:128),I3(1:128,129:256);I3max()(129:256,1:128),I(129:256,129:256)];%25%shi,threshold is 0.5263,output is 0.5544
%figure,imshow(I3)
%-------------------------------------------------
%J2=imcrop(I3,[64 64 127 127]);%1/4
%J2=imcrop(I3,[64 64 63 63]);%1/8
%J3=[zeros(128,63),J2,zeros(128,65)];%1/4
%J3=[zeros(64,63),J2,zeros(64,129)];%1/8
%I3=[zeros(63,256);J3;zeros(65,256)];%1/4
%I3=[zeros(63,256);J3;zeros(129,256)];%1/8
%------------------------------------------------
%resizing
%------------------------------------------------
J=imresize(I3,.5,'nearest');
I3=imresize(J,2,'nearest');
%------------------------------------------------
I3=double(I3)/255;
%figure,imshow(I3,[])
invariancematrix3= abs(invariance(I3*255,1,0));   %原始圖像的不變矩
%------------------------------------------------
MSE=sqrt(sum(sum((I3*255-I*255).^2))/prod(size(I)));%均方根誤差(I'-I)^2/N^2;
PSNR=20*log10(255/MSE);%水印算法中的峰峰信噪比，PSNR=10log(255^2/MSE^2),35.5125
PSNR
%------------------------------------------------
Di=blkproc(I3,[8 8],'P1*x*P2',dctm,dctm');
Dd=blkproc(Di,[8 8],'chaoswm(x)');
%Mud=mean2(Dd)
%STDd=std2(Dd)
%Dd=(Dd-min(min(Dd)))/(max(max(Dd))-min(min(Dd)));
%------------------------------------------------
x1=linspace(0.2000000000000,0.2000000000002,1001);%
for i=1:1001,
   w1=choasq(x1(i),128);
   w2=sum(sum(Dd.^2));
   w3=sum(sum((w1).^2));
   Lou(i)=sum(sum(Dd.*w1))/w2;
end
%------------------------------------------------
%w1=msqu(127,7);
%w2=sum(sum(Dd.^2));
%w3=sum(sum(w1.^2));
%Lou=sum(sum(Dd.*w1))/w2
%------------------------------------------------
Mul=mean2(Lou);
STDl=std2(Lou);
Lou=Lou-Mul;
max(Lou);
%Thres=L(999);
x2=1:1000;
y1=ones(1,1000).*(4*STDl);%自動閾值的選取，'4σ'原則
y2=ones(1,1000).*(-4*STDl);
%y1(1),y2(1)
figure,plot(x2,Lou(2:1001),'-',x2,y1,':',x2,y2,':');
xlabel('Watermarks')
ylabel('Detector Response')
legend('responds','threshold')
invariantmatrx=[invariancematrix1;invariancematrix2;invariancematrix3]
y1=sum((invariancematrix1-invariancematrix3).^2)