数字水印matlab程序.docx
《数字水印matlab程序.docx》由会员分享，可在线阅读，更多相关《数字水印matlab程序.docx（27页珍藏版）》请在冰豆网上搜索。
数字水印matlab程序
clearall;
closeall;
clc;
V=double（imread（'D:
\lena\lena.jpg'））;
imshow（mat2gray（V））;
[iu]=size（V）;%计算V的规格
r=100;%设置分解矩阵的秩
W=rand（i,r）%初始化WH，为非负数
H=rand（r,u）
maviter=100;%最大迭代次数
foriter=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:
Gaussiannoise'）;
subplot（2,2,3）;
title（'edgeoforigineimage'）
[Y0,I0]=sort（watermark1）;
%
%读入原图象
trueImage=imread（'H:
\DocumentsandSettings\sunhw\MyDocuments\MyPictures\biaozhun.bmp'）;
alfa=.1;
LENGTH=2500;
subplot（2,2,2）;
imshow（trueImage）;
title（'origineimage:
I'）;
%
%对原图象进行DCT变换
dctF1=dct2（'H:
\DocumentsandSettings\sunhw\MyDocuments\MyPictures\biaozhun.bmp'）;
[m,n]=size（dctF1）;
%
%找出水印嵌入位置（幅值较大的n个频域成分）
A=dctF1（:
）;
[Y1,I1]=sort（A）;
x=m*n;
k=LENGTH;
M=zeros（x,1）;
%
%修改幅值较大的n个频域成分的幅值，嵌入水印（因为两个问题不同，所以有两个注释符）
fori=1:
x
ifk>=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;
fori=1:
x
N（I1（i））=M（i）;
end
a=1;
forj=1:
n
fori=1:
m
dctF2（i,j）=N（a）;
a=a+1;
end
end
%
%DCT反变换，得到嵌入水印的图象
idctF1=idct2（dctF2）;
subplot（2,2,4）;
imshow（log（abs（idctF1））,[]）;
title（'embededimage:
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'）;
[iu]=size（V）;%计算V的规格
r=100;%设置分解矩阵的秩
W=rand（i,r）%初始化WH，为非负数
H=rand（r,u）
maviter=100;%最大迭代次数
foriter=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
fori=1:
rm %i=1:
32
forj=1:
cm %j=1:
32
x=（i-1）*10;y=（j-1）*10;
ifmark（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
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）;
form=1:
InfoStrSize
Infochar=double（Info（m））; %%'c'为99
forn=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;
form=1:
blockno
forn=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）;%得到边界矩阵。
ifm<=1|n<=1
T=0;
else
T=sum（BW_8_8）; T=sum（T）;
end
ifT>T1
Alpha=Alpha2;
%block_dct1（1,1）=block_dct1（1,1）*（1+Alpha*mark（k））;
ifl<=（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
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;
form=1:
blockno
forn=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）;%得到边界矩阵。
ifm<=1|n<=1
T=0;
else
T=sum（BW_8_8）; T=sum（T）;
end
ifT>T1
Alpha=Alpha2;
%block_dct1（1,1）=block_dct1（1,1）*（1+Alpha*mark（k））;
ifl<=（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
end
subplot（3,2,5）;
plot（array2）;
title（'提取水印'）;
extractedInfo=zeros（InfoStrSize,1）;
form=1:
InfoStrSize
infochar=0;
forn=1:
8
ifarray2（8*（m-1）+n）==1
infochar=infochar+bitset（0,n,1）;
end
end
extractedInfo（m）=infochar+extractedInfo（m）;
end
resultStr=char（extractedInfo）;
subplot（3,2,6）;
plot（array2）;
r=Arnold（w0,row,colum,times）
fork=1:
times
fori=1:
row
forj=1:
colum
i1=i+j;
j1=i+2*j;
ifi1>row
i1=mod（i1,row）;
end
ifj1>colum
j1=mod（j1,colum）;
end
ifi1==0
i1=row;
end
ifj1==0
j1=colum;
end
w1（i1,j1）=w0（i,j）;
end
end
w0=w1;
end
r=w0;
%提取水印
forp=1:
size/B
forq=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
end
end
figure（5）;imshow（F,[]）;title（'提取出的水印'）;
%
%攻击实验
disp（'inputyouchoiceaccordingtothefollowing
imageprocessingoperation:
'）;
disp（'0--exit'）;
disp（'1--smoothingpatterns'）;
%添加噪音
disp（'2--addinguniormnoise添加噪音'）;
%滤波
disp（'3--addingfilter[1010]滤波'）;
%剪切
disp（'4--cuttingpartoftheimage剪切'）;
%压缩
disp（'5--10qualityJPEGcompressing压缩'）;
%旋转45度
disp（'6--rotate45旋转'）;
%
d=input（'pleaseinputyouchoice（1,2,3,4,5,6）:
'）;
whiled~=0
switchd
case1
watermark_detect（idctF1,Y1,I0,waterMark1）;
case2
WImage2=idctF1;
noise0=10*rand（size（WImage2））;
WImage2=WImage2+noise0;
figure;
imshow（WImage2,[]）;
title（'addinguniformnoise添加噪音'）;
watemark_detect（WImage2,Y1,I0,waterMark1）;
case3
WImage3=idctF1;
H=fspcial（'gaussian高斯',[10,10],5）;
WImage3=imfilter（WImage3,H）;
figure;
imshow（WImage3,[]）;
title（throughfilter[10,10]滤波'）;
watemark_detect（WImage3,Y1,I0,waterMark1）;
case4
WImage4=idctF1;WImage4（1:
128,1;128）=256;
figure;
imshow（WImage4）;
title（'cuttingpartoftheimage剪切'）;
watemark_detect（WImage4,Y1,I0,waterMark1）;
case5
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（'JPEGImage压缩'）;
watemark_detect（WImage5,Y1,I0,waterMark1）;
case6WImage6=idctF1;
WImage6=imrotate（WImage6,45,'bilinear','corp'）;
figure;
imshow（Wimage6）;
title（'rotate45旋转'）;
watemark_detect（WImage6,Y1,I0,waterMark1）;
case0
break;
otherwise
error（'youhaveavalidvalue（您的输入错误）'）;
end
d=input（'pleaseinputyouchoice（请输入您的选择）:
'）;
end
%结束
f=imread（'watermark.jpg'）;
%将含水印图像f归一化,以便于攻击处理。
m=max（max（I_W）;
I_W=double（I_W）./double（m）;
%---攻击-------------------------------------------------------------------
attack=0;
switchattack
case0,
attackf=I_W;
att='未攻击';
case1,
%%1.JPEG压缩
imwrite（I_W,'attackf.jpg','jpg','quality',30）;
attackf=imread（'attackf.jpg'）;
attackf=double（attackf）/255;
att='JPEG压缩';
case2,
%%2.高斯低通滤波
h=fspecial（'gaussian',3,1）;
attackf=filter2（h,I_W）;
att='高斯低通滤波';
case3,
%%3.直方图均衡化
attackf=histeq（I_W）;
att='直方图均衡化';
case4,
%%4.图像增亮
attackf=imadjust（I_W,[],[0.4,1]）;
att='图像增亮';
case5,
%%5.图像变暗
attackf=imadjust（I_W,[],[0,0.85]）;
att='图像变暗';
case6,
%%6.增加对比度
attackf=imadjust（I_W,[0.3,0.6],[]）;
att='增加对比度';
case7,
%%7.降低对比度
attackf=imadjust（I_W,[],[0.2,0.8]）;
att='降低对比度';
case8,
%%8.添加高斯噪声
attackf=imnoise（I_W,'gaussian',0,0.01）;
att='添加高斯噪声';
case9,
%%9.椒盐噪声
attackf=imnoise（I_W,'salt&pepper',0.06）;
att='椒盐噪声';
case10,
%%10.添加乘积性噪声
attackf=imnoise（I_W,'speckle',0.08）;
att='添加乘积性噪声';
end;
%---攻击后处理--------------------------------------------------------------
I_W=attackf.*double（m）;
figure
（2）;
imshow（uint8（I_W）;%显示水印嵌入图攻
升级会员 