当前位置：文档之家› 数字水印matlab程序

数字水印matlab程序

clear all;

close all;

clc;

V=double(imread('D:\lena\lena.jpg'));

imshow(mat2gray(V));

[i u]=size(V); %计算V的规格

r=100; %设置分解矩阵的秩

W=rand(i,r) %初始化WH，为非负数

H=rand(r,u)

maviter=100; %最大迭代次数

for iter=1:maviter

W=W.*((V./(W*H))*H'); %注意这里的三个公式和文中的是对应的

W=W./(ones(i,1)*sum(W));

H=H.*(W'*(V./(W*H)));

end

img_V=W*H;

figure;

imshow(mat2gray(img_V));

首先读入原始图象并设置参数，然后嵌入水印信息，程序代码如下：clear

size=256;block=8;blockno=size/block;LENGTH=size*size/64;

Alpha1=0.02; Alpha2=0.1; T1=3; I=zeros(size,size); D=zeros(size,size);

BW=zeros(size,size); block_dct1=zeros(block,block);

%产生水印序列并对其排序

randn('seed',10);watermark1=randn(1,LENGTH);

subplot(2,2,1);plot(watermark1);title('watermarc:Gaussian noise');

subplot(2,2,3);

title('edge of origine image')

[Y0,I0]=sort(watermark1);

%读入原图象

trueImage=imread('H:\Documents and Settings\sunhw\My Documents\My Pictures\biaozhun.bmp');

alfa=.1;

LENGTH=2500;

subplot(2,2,2);

imshow(trueImage);

title('origine image:I');

%对原图象进行DCT变换

dctF1=dct2('H:\Documents and Settings\sunhw\My Documents\My

Pictures\biaozhun.bmp');

[m,n]=size(dctF1);

%找出水印嵌入位置（幅值较大的n个频域成分）

A=dctF1(:);

[Y1,I1]=sort(A);

x=m*n;

k=LENGTH;

M=zeros(x,1);

%修改幅值较大的n个频域成分的幅值，嵌入水印（因为两个问题不同，所以有两个注释符）

for i=1:x

if k>=1

M(x)=Y1(x)*(1+alfa*Y0(k));

k=k-1;

else

M(x)=Y1(x);

end

x=x-1;

end

N=zeros(x,1);

x=m*n;

for i=1:x

N(I1(i))=M(i);

end

a=1;

for j=1:n

for i=1:m

dctF2(i,j)=N(a);

a=a+1;

end

%DCT反变换，得到嵌入水印的图象

idctF1=idct2(dctF2);

subplot(2,2,4);

imshow(log(abs(idctF1)),[ ]);

title('embeded image:D');

%end

I=imread('D:\lena\1.jpg');

disp(I);

I=double(I)/255;

disp(I);

I=ceil(I);

%%%%%%%%%%显示水印图像%%%%%%%%%%%%%

figure(1);

subplot(2,3,1);

imshow(I),title('水印图像')

dimI=size(I);

rm=dimI(1);cm=dimI(2);

%%%%%%%%%%%%%%%5 以下生成水印信息 %%

mark=I;

alpha=0.05;

V=imread('D:\lena\lena.jpg');

[i u]=size(V); %计算V的规格

r=100; %设置分解矩阵的秩

W=rand(i,r) %初始化WH，为非负数

H=rand(r,u)

maviter=100; %最大迭代次数

for iter=1:maviter

W=W.*((V./(W*H))*H'); %注意这里的三个公式和文中的是对应的 W=W./(ones(i,1)*sum(W));

H=H.*(W'*(V./(W*H)));

end

k1=H;

psnr_cover=double(V);

subplot(2,3,2),imshow(a0,[]),title('载体图像');

[r,c]=size(a0);

cda0=blkproc(a0,[8,8],'dct2'); %%%%%%%%%%%%%%%%%%%% 嵌入 %%%%%% cda1=cda0; % cda1 = 256_256

for i=1:rm % i=1:32

for j=1:cm % j=1:32

x=(i-1)*10;y=(j-1)*10;

if mark(i,j)==1

k=k1;

else

k=k2;

end

cda1(x+1,y+8)=cda0(x+1,y+8)*(1+alpha*k(1));

cda1(x+2,y+7)=cda0(x+2,y+7)*(1+alpha*k(2));

cda1(x+3,y+6)=cda0(x+3,y+6)*(1+alpha*k(3));

cda1(x+4,y+5)=cda0(x+4,y+5)*(1+alpha*k(4));

cda1(x+5,y+4)=cda0(x+5,y+4)*(1+alpha*k(5));

cda1(x+6,y+3)=cda0(x+6,y+3)*(1+alpha*k(6));

cda1(x+7,y+2)=cda0(x+7,y+2)*(1+alpha*k(7));

cda1(x+8,y+1)=cda0(x+8,y+1)*(1+alpha*k(8));

end

%%%%% 嵌入水印后图像 %%%%%%%%%%%%%%

a1=blkproc(cda1,[8,8],'idct2');

a_1=uint8(a1);

imwrite(a_1,'withmark.bmp','bmp');

subplot(2,3,3),imshow(a1,[]),title('嵌入水印后的图像');

size=256;

block=8;

blockno=size/block;%一行有32格

LENGTH=size*size/64;

Alpha1=0.025;

Alpha2=0.1;

T1=3;

I=zeros(size,size);%产生全矩阵

D=zeros(size,size);

BW=zeros(size,size);

Block_dct1=zeros(block,block);

%产生水印，并显示水印信息；

subplot(3,2,1);

Info='dcf';

InfoStrSize=length(Info);

%将字符串转换为位数组

array=zeros(1,InfoStrSize*8);

for m=1:InfoStrSize

Infochar=double(Info(m)); %% 'c'为99

for n=1:8

array(8*(m-1)+n)=bitget(Infochar,n);%%获得Infochar第n位的值 end

end

plot(array);

title('原始水印信息');

%显示原图

subplot(3,2,2);

i=imread('lena.bmp');

imshow(i,[]);

title('原始图像')

%显示prewitt为算子的边缘图

BW=edge(i,'prewitt');

%BW=edge(I,’Roberts’);

%BW=edge(I,’Sobel’);

%BW=edge(I,’zerocross’);

subplot(3,2,3);imshow(BW);

Title('原始图像边缘图');

%嵌入水印

l=1;

k=1;

for m=1:blockno

for n=1:blockno

x=(m-1)*block+1; y=(n-1)*block+1;%算出每格图像的坐标（x,y）,block=8,8*8的图像小格

block_dct1=H(x:x+block-1,y:y+block-1);%取原始图像小格中的像素点到block_dct1矩阵中。

block_dct1=dct2(block_dct1);%对二维数组进行离散余弦变换。dct是有损压缩如jpeg使用的技术。Dct是可逆的运算

BW_8_8=BW(x:x+block-1,y:y+block-1);%得到边界矩阵。

if m<=1|n<=1

T=0;

else

T=sum(BW_8_8); T=sum(T);

end

if T>T1

Alpha=Alpha2;

%block_dct1(1,1)=block_dct1(1,1)*(1+Alpha*mark(k));

if l<=(InfoStrSize*8)

block_dct1(1,1)=block_dct1(1,1)*(1+Alpha*array(l));

l=l+1;

end

else

Alpha=Alpha1;

% block_dct1(1,1)=block_dct1(1,1)*(1+Alpha*mark(k));

end

Block_dct=idct2(block_dct1);

D(x:x+block-1,y:y+block-1)=Block_dct;

k=k+1;

end

%显示嵌入水印后的图像

subplot(3,2,4);imshow(D,[]);title('嵌入水印图像')

%保存该图像

D=uint8(D);

imwrite(D,'marked.bmp');

2：

%提取水印

D=imread('marked.bmp');

D=double(D);

I=imread('lena.bmp');

I=double(I);

array2=zeros(1,InfoStrSize*8);

K=1;

l=1;

for m=1:blockno

for n=1:blockno

x=(m-1)*block+1; y=(n-1)*block+1;%算出每格图像的坐标（x,y）,block=8,8*8的图像小格

block_dct1=I(x:x+block-1,y:y+block-1);%取原始图像小格中的像素点到block_dct1矩阵中。

block_dct2=D(x:x+block-1,y:y+block-1);

Block_dct1=dct2(block_dct1);%对二维数组进行离散余弦变换。dct是有损压缩如jpeg使用的技术。Dct是可逆的运算

Block_dct2=dct2(block_dct2);

BW_8_8=BW(x:x+block-1,y:y+block-1);%得到边界矩阵。

if m<=1|n<=1

T=0;

else

T=sum(BW_8_8); T=sum(T);

end

if T>T1

Alpha=Alpha2;

%block_dct1(1,1)=block_dct1(1,1)*(1+Alpha*mark(k));

if l<=(InfoStrSize*8)

tmp=(Block_dct2(1,1)/Block_dct1(1,1)-1);

tmp=tmp/Alpha;

tmp2=round(tmp);

array2(l)=double(tmp2);

l=l+1;

end

else

Alpha= Alpha1;

%block_dct1(1,1)=block_dct1(1,1)*(1+Alpha*mark(k));

end

k=k+1;

end

subplot(3,2,5);

plot(array2);

title('提取水印');

extractedInfo=zeros(InfoStrSize,1);

for m=1:InfoStrSize

infochar=0;

for n=1:8

if array2(8*(m-1)+n)==1

infochar=infochar+bitset(0,n,1);

end

extractedInfo(m)=infochar+extractedInfo(m);

end

resultStr=char(extractedInfo);

subplot(3,2,6);

plot(array2);

r=Arnold(w0,row,colum,times)

for k=1:times

for i=1:row

for j=1:colum

i1=i+j;

j1=i+2*j;

if i1>row

i1=mod(i1,row);

end

if j1>colum

j1=mod(j1,colum);

end

if i1= =0

i1=row;

end

if j1= =0

j1=colum;

end

w1(i1,j1)=w0(i,j);

end

w0=w1;

end

r=w0;

%提取水印

for p=1:size/B

for q=1:size/B

x=(p-1)*B+1;

y=(q-1)*B+1;

if (I_W(x,y)-P(x,y))>0

F(p,q)=1;

else

F(p,q)=0;

end

figure(5); imshow(F,[]); title('提取出的水印');

%攻击实验

disp('input you choice according to the following

image processing operation:');

disp('0--exit');

disp('1--smoothing patterns');

%添加噪音

disp('2--adding uniorm noise 添加噪音');

%滤波

disp('3--adding filter [10 10] 滤波');

%剪切

disp('4--cutting part of the image 剪切');

%压缩

disp('5--10 quality JPEG compressing 压缩');

%旋转45度

disp('6--rotate 45 旋转');

d=input('please input you choice(1,2,3,4,5,6):');

while d~=0

switch d

case 1

watermark_detect(idctF1,Y1,I0,waterMark1);

case 2

WImage2=idctF1;

noise0=10*rand(size(WImage2));

WImage2=WImage2+noise0;

figure;

imshow(WImage2,[ ]);

title('adding uniform noise 添加噪音');

watemark_detect(WImage2,Y1,I0,waterMark1);

case 3

WImage3=idctF1;

H=fspcial('gaussian高斯',[10,10],5);

WImage3=imfilter(WImage3,H);

figure;

imshow(WImage3,[ ]);

title(through filter [10,10] 滤波');

watemark_detect(WImage3,Y1,I0,waterMark1);

case 4

WImage4=idctF1; WImage4(1:128,1;128)=256; figure;

imshow(WImage4);

title('cutting part of the image 剪切');

watemark_detect(WImage4,Y1,I0,waterMark1);

case 5

WImage5=idctF1;

WImage5=im2double(WImage5);

cnum=10;

dctm=dctmtx(8);

p1=dctm;

p2=dctm.';

imageDCT=blkproc(WImage5,[8,8],'p1*p2*x',dctm,dctm.');

DCTvar=im2col(imageDCT,[8,8],'distinct').';

n=size(DCTvar,1);

DCTvar=(sum(DCTvar.*DCTvar)-(sum(DCTvar)/n).^2)/n;

[dum,order]=sort(DCTvar);

cnum=64-cnum;

mask=ones(8,8);

mask(order(1:cnum))=zeros(1,cnum);

im88=zeros(9,9);

im88(1:8,1:8)=mask;

im128128=kron(im88(1:8,1:8),ones(16));

dctm=dctmtx(8);

p1=dctm.';

p2=mask(1;8,1:8);

p3=dctm;

Wimage5=bikproc(imageDCT,[8,8],'p1*(x.8p2)*p3',dctm.',mask(1:8,1:8),dctm); figure;

imshow(Wimage5);

title('JPEG Image 压缩');

watemark_detect(WImage5,Y1,I0,waterMark1);

case 6 WImage6=idctF1;

WImage6=imrotate(WImage6,45,'bilinear','corp');

figure;

imshow(Wimage6);

title('rotate 45 旋转');

watemark_detect(WImage6,Y1,I0,waterMark1);

case 0

break;

otherwise

error('you have a valid value(您的输入错误)');

end

d=input('please input you choice(请输入您的选择):');

end

%结束

f=imread('watermark.jpg');

%将含水印图像f归一化,以便于攻击处理。

m=max(max(I_W);

I_W=double(I_W )./double(m);

%---攻击-------------------------------------------------------------------

attack=0;

switch attack

case 0,

attackf= I_W;

att='未攻击';

case 1,

%%1. JPEG 压缩

imwrite(I_W,'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, I_W);

att='高斯低通滤波';

case 3,

%%3. 直方图均衡化

attackf=histeq(I_W);

att='直方图均衡化';

case 4,

%%4. 图像增亮

attackf=imadjust(I_W,[],[0.4,1]);

att='图像增亮';

case 5,

%%5. 图像变暗

attackf=imadjust(I_W,[],[0,0.85]);

att='图像变暗';

case 6,

%%6. 增加对比度

attackf=imadjust(I_W,[0.3,0.6],[]);

att='增加对比度';

case 7,

%%7. 降低对比度

attackf=imadjust(I_W,[],[0.2,0.8]);

att='降低对比度';

case 8,

%%8. 添加高斯噪声

attackf=imnoise(I_W,'gaussian',0,0.01); att='添加高斯噪声';

case 9,

%%9. 椒盐噪声

attackf=imnoise(I_W,'salt & pepper',0.06); att='椒盐噪声';

case 10,

%%10. 添加乘积性噪声

attackf=imnoise(I_W,'speckle',0.08);

att='添加乘积性噪声';

end;

%---攻击后处理--------------------------------------------------------------

I_W=attackf.*double(m);

figure(2);

imshow(uint8(I_W);%显示水印嵌入图攻击后效果

title(att);

imwrite(uint8(I_W),'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;

V=imread('D:\lena\lena.jpg');

Alfa=0.05;

LENGTH=2500;

[i u]=size(V); %计算V的规格

r=100; %设置分解矩阵的秩

W=rand(i,r) %初始化WH，为非负数

H=rand(r,u)

maviter=100; %最大迭代次数

for iter=1:maviter

W=W.*((V./(W*H))*H'); %注意这里的三个公式和文中的是对应的

W=W./(ones(i,1)*sum(W));

H=H.*(W'*(V./(W*H)));

end

[Y1,10]=sort(H)

L=imread('D:\lena\1.jpg');

[Y0,10]=sort(L)

x=m*n;

k=LENGTH;

M=zeros(x,1);

%修改幅值较大的n个频域成分的幅值，嵌入水印（因为两个问题不同，所以有两个

注释符）

for i=1:x

if k>=1

M(x)=Y1(x)*(1+alfa*Y0(k));

k=k-1;

else

M(x)=Y1(x);

end

x=x-1;

end

N=zeros(x,1);

x=m*n;

for i=1:x

N(I1(i))=M(i);

end

a=1;

for j=1:n

for i=1:m

N(i,j)=N(a);

a=a+1;

end

%重构图像，得到嵌入水印的图象B=W.*H;

figure(1);

imshow(B);

title('嵌入水印后的图象');

end

clear all;

clc;

start_time=cputime; %%%%%%%%%%%% 读取水印图像 %%%%%%%% I=imread('mark.bmp');

I=rgb2gray(I);

I=double(I)/255;

I=ceil(I);

%%%%%%%%%%显示水印图像%%%%%%%%%%%%% figure(1);

subplot(2,3,1);

imshow(I),title('水印图像')

dimI=size(I);

rm=dimI(1);cm=dimI(2);

%%%%%%%%%%%%%%%5 以下生成水印信息 %%

mark=I;

alpha=50,

k1=randn(1,8);

k2=randn(1,8);

a0=imread('lena.bmp');

psnr_cover=double(a0);

subplot(2,3,2),imshow(a0,[]),title('载体图像');

[r,c]=size(a0);

cda0=blkproc(a0,[8,8],'dct2'); %%%%%%%%%%%%%%%%%%%%% 嵌入 %%%%%%%%%%

cda1=cda0; % cda1 = 256_256

for i=1:rm % i=1:32

for j=1:cm % j=1:32

x=(i-1)*8;y=(j-1)*8;

if mark(i,j)==1

k=k1;

else

k=k2;

end

cda1(x+1,y+8)=cda0(x+1,y+8)+alpha*k(1);

cda1(x+2,y+7)=cda0(x+2,y+7)+alpha*k(2);

cda1(x+3,y+6)=cda0(x+3,y+6)+alpha*k(3);

cda1(x+4,y+5)=cda0(x+4,y+5)+alpha*k(4);

cda1(x+5,y+4)=cda0(x+5,y+4)+alpha*k(5);

cda1(x+6,y+3)=cda0(x+6,y+3)+alpha*k(6);

cda1(x+7,y+2)=cda0(x+7,y+2)+alpha*k(7);

cda1(x+8,y+1)=cda0(x+8,y+1)+alpha*k(8);

end

%%%%% 嵌入水印后图像 %%%%%%%%%%%%%%

a1=blkproc(cda1,[8,8],'idct2');

a_1=uint8(a1);

imwrite(a_1,'withmark.bmp','bmp');

subplot(2,3,3),imshow(a1,[]),title('嵌入水印后的图像');

disp('嵌入水印处理时间');

embed_time=cputime-start_time,

%%%%%%%%

disp('对嵌入水印的图像的攻击实验，请输入选择项：'); disp('1--添加白噪声');

disp('2--高斯低通滤波');

disp('3--JPEG 压缩');

disp('4--图像剪切');

disp('5--旋转10度');

disp('6--直接检测水印');

disp('其他--不攻击');

d=input('请输入选择（1-6）:');

start_time=cputime;

figure(1);

switch d

case 6

subplot(2,3,4);

imshow(a1,[]);

title('未受攻击的含水印图像');

M1=a1;

case 1

WImage2=a1;

noise0=20*randn(size(WImage2));

WImage2=WImage2+noise0;

subplot(2,3,4);

imshow(WImage2,[]);

title('加入白噪声后图像');

M1=WImage2;

M_1=uint8(M1);

imwrite(M_1,'whitenoise.bmp','bmp');

case 2

WImage3=a1;

H=fspecial('gaussian',[4,4],0.2);

WImage3=imfilter(WImage3,H);

subplot(2,3,4);

imshow(WImage3,[]);