1、第三章 图像数字水印的方案3.1 图像数字水印的技术方案在数据库中存储在国际互联网上传输的水印图像一般会被压缩，有时达到很高的压缩比。因此，数字水印算法所面临的第一个考验就是压缩。JPEG 和EZW（Embedded Zero-Tree Wavelet）压缩是最常见的两种压缩方法。 JPEG 是基于离散余弦变换域的压缩方法，而 EZW 是基于小波变换域的压缩方法。前人的研究证明采用与压缩算法相同的变换域水印方法，对于压缩的稳健性较强。因此，我研究图像文件水印算法主要集中在变换域算法及利用人眼视觉特性上。数字水印的嵌入要求即要考虑视觉透明性，又要保证嵌入水印后图像的稳健性，这两个方面存在着矛盾。
2、保证视觉透明性，就要将水印嵌入到人眼不敏感区，也就是嵌入到图像的高频分量中。而多数图像处理方法对于图像高频部分的损坏程度较高，如有损压缩、高频滤波等。水印很容易在经历图像处理的过程中丢失。这样，则无法保证图像数字水印的稳健性。如果要获得很好的稳健性，数字水印应加在人眼敏感的低频部分，图像的大部分能量集中在低频部分，如果对于低频部分进行处理，水印固然会失去，而图像也没有了利用价值，然而，水印的嵌入会对图像的质量有非常大的影响，这又无法保证视觉透明性。数字水印算法的实现基本分为三个部分：宿主图像的变换，水印的嵌入和水印的检测，分别描述如下。3.2 基于 DCT 域的图像数字水印技术离散余弦变换（D
3、iscrete Cosine Transform）属于正交变换图像编码方法中的一种。正交变换图像编码始于 1968 年。当时安德鲁斯（Andrews）等人发现大多数自然图像的高频分量相对幅度较低，可完全舍弃或者只用少数码字编码，提出不对图像本身编码，只对其二维傅立叶（DFT）系数进行编码和传输。但 DFT 是一种正交变换，运算量很大，常常使实时处理发生困难，第二年他们就用 Walsh-Hadamard 变换（WHT）取代 DFT 可以使运算量明显减少，这是因为 WHT 变换只有加减法而无需乘法。但是更有意义的是离散余弦变换和离散正旋变换的出现，它们具有快速算法，精确度高。其中最重要的是 197
4、4 年提出的 DCT，因为其变换矩阵的基向量很近似于托伯利兹矩阵的特征向量，而托伯利兹矩阵又体现了人类语言及图像信号的相关性。因此，DCT 常常被认为是语音与图像信号变换的准最佳变换。图像是二维的，所以在研究时主要用到二维 DCT，以及二维 IDCT 来对图像进行处理。3.2.1 离散余弦变换（DCT）的定义数字图像 X（m,n）是具有 M 行 N 列的一个矩阵。为了同时减弱或去除图像数据相关性，可以运用二维 DCT，将图像从空间域转换到 DCT 变换域。根据定义，二维离散余弦变换（DCT）定义如下： NlnkmnXlckNlkYMmNn 2)1(cos2)1(cos),()(2),(10 式
5、中：m，k=0,1,M-1; n， =0,1,N-1。l1,210)( Mkkc其 中 函 数,210)( Nlc二维逆离散余弦变换（IDCT）的定义如下： 10 2)1(cos2)1(cos),(),(MKNL NlnMkmlkYcnmX 式中：m，k=0,1,M-1; n， =0,1,N-1。l3.2.2 离散余弦变换的特点在基于 DCT 的变换编码中，图像是先经分块（8 8 或 1616）后再经 DCT，这种变换是局部的，只反映了图像某一部分的信息。当然也可以对整幅图像的特点，但是运算速度比分块 DCT 要慢。图像经 DCT 后，得到的 DCT 图像有三个特点：一是系数值全部集中到 0
6、值附近（从直方图统计的意义上） ，动态范围很小，这说明用较小的量化比特数即可表示 DCT 系数；二是 DCT 变换后图像能量集中在图像的低频部分，即 DCT 图像中不为零的系数大部分集中在一起（左上角） ，因此编码效率很高。三是没有保留原图像块的精细结构，从中反映不了原图像块的边缘、轮廓等信息，这一特点是由 DCT 缺乏时局域性造成的。如下左图 31 是原始图像经过 DCT 变换后的系数图像为图 32。两条线划分出图像的低频、中频和高频分别所在的矩形区域。可以看出，图像 DCT 变换后大部分参数接近于零，只有左上角的低频部分有较大的数值，中频部分参数值相对较小，而大部分高频参数值非常小，接近于
7、零。图 31：原图像 图 32：变换后的系数图像3.2.3 离散余弦变换的数字水印算法根据离散余弦变换后的参数性质，本文采用了以 ZigZag 方式重排变换域系数的方法，选出中频分量，用数字水印序列对其进行非线性调制。水印检测时，待检测图像仍按比方式选择变换域系数，与待水印进行相关运算，与阈值比较来判断是否所含水印。离散余弦域的数字水印算法的具体实现分为三步：宿主图像的变换，数字水印的嵌入，数字水印的检测。3.2.3.1 宿主图像的 DCT 变换对于 NN 大小的 256 灰度级的宿主图像 I 进行 NN 二维离散余弦变换（ DCT） 。以 ZigZag 方式对于 DCT 变换后的图像频率系数
8、重新排列成一维向量 Y=y1, y2,yNN.并取出序列中第 L+1 到 L+M 的中频系数部分,得到 YL= YL+1, YL+2, YL+M 3.2.3.2 数字水印的嵌入假设数字水印 W 为一服从标准正态分布的随机实数序列，用数字序列表示为W=W1 ，W 2 ，W M 。用 W 对 Y 序列中第 L+1 到 L+M 的中频系数部分的值进行修改，按以下公式进行： 2,0iiiyiLmw或经过修改的系数序列 Y = Y 1, Y 2,. Y NN 以 ZigZag 逆变换形式重组，再进行 NN DCT 逆变换，得到嵌有数字水印的图像 I。 .3.2.3.3 数字水印的检测待检测的可能含有水印
9、的图像 I“ 。假设 I“未损失大量信息，可以近似认为 I“= I。在此假设下可以运用统计的方法来检测水印。（1 ）待检水印域待检图像中频系数相关性的测定同样对 I进行 DCT 变换，以 ZigZag 方式将 DCT 系数排成一维向量 Y “= Y1 “, Y2 “,. YNN“。由于假设 I“=I，则 Y“= Y 。取出 Y“（等于 Y）中第 L+1 到 L+M 的中频系数部分 YL “= YL+1 “, YL+2“,. YL+M“。假设待检测的数字水印 X=X1，X 2，. XM为一符合标准正态分布的实数伪随机序列。则可以通过待检水印与图像中频系数作相关运算来判断是否所加入了水印。只有在待
10、检水印为所加入的水印时，才能得到较大的相关值。否则相关值很小，接近于零。 )()(1)(1),cov( 2XiWYXYMXYZ iLiLiiLL 用符号 E 表示数学期望，得到：2()()0zYWX没 有 水 印 存 在（2）阈值的确定根据中心极限定理，参照水印匹配与不匹配两种情况得到阈值为 。2/)(LzYET由于原始图像难以得到，因此从实用性出发，阈值定义为： MiiLLZYYET12/)(综上所述，满足 时，则表明检测到匹配水印。否则，未检测5.02ZZT或到匹配的水印。3.3 MATLAB 工具简介3.3.1. 简介Matlab 是当前在国内外十分流行的工程设计和系统仿真软件包。它是
11、MathWorks公司于 1982 年推出的一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算、信号处理和图形显示于一体，构成了一人方便的、界面友好的用户环境。Matlab 的推出得到了各个领域专家、学者的广泛关注，其强大的扩展功能为各个领域的应用提供了基础。由各个专家学者相继推出了 MATLAB 工具箱，其中的信号处理(signal processing)、控制系统(control system)、神经网络 (neural network)、图像处理(image processing)、鲁棒控制(robust control)、非线性系统控制设计(nonlinear system co
12、ntrol design)、系统辨识(system identification)、最优化(optimization)、模糊逻辑(fuzzy logic)、小波(wavelet)、通信(communication)、统计(statistics)等工具箱，这些工具箱给各个领域的研究和工程应用提供了有力的工具，借助于这些“巨人肩上的工具” ，各个层次的研究人员可直观、方便地进行分析、计算及设计工作，从而大大地节省了时间。3.3.2. MATLAB 研究数字水印的优点 集成了 DCT、DWT 等函数有丰富的小波函数和处理函数，这不仅方便了研究人员，而且使源程序简洁明了、易实现。 强大的数学运算功能。
13、能够方便、高效地实现音频、视频中的大量矩阵运算。 提供了图像处理工具箱、小波分析工具箱、数字信号处理工具箱。用来编制跨数字图像处理技术、数字信号处理等多学科的数字水印技术是非常好的选择。 MATLAB 与目前最强大的编程工具 Visual C+具有良好的接口。3.3.3. MATLAB 函数介绍在介绍函数之前，我们必须明确一点：作水印程序时，处理的图像数据是二维信号，而声音信号是一维信号。这里，我们仅仅简单介绍与水印有关的函数。 数据输入输出函数imread（）和 imwrite（）：可以读写 bmp,jpg/jpeg, tif/tiff, png, hdf, pcx, wxd 格式文件。读索
14、引文件时，还可以得到相应的调色板数据。auread（） 、auwrite（） 、wavread（）和 wavwrite（）：可以方便地读写 au 和 wav 文件，并可控制其中的位及频率。 图像显示imshow()：显示一幅图像；imfinfo()：可以得到读入图像的信息。如文件的大小、格式、格式版本号、图像的高度、宽度、颜色类型（真彩色，灰度图还是索引图）等。 变换频函数对信号采用不同的变换，是实现频域法水印的至关重要的一步，MATLAB 中对一维信号和二维信号分别提供了各种变换和逆变换函数。 离散余弦变换（DCT）dct(),dct2():分别实现一维信号和二维信号的 DCT（离散余弦变换
15、） ；idct(),idct2():分别实现一维信号和二维信号的 IDCT（逆向离散余弦变换） ； 离散小波变换（DWT）dwt(),dwt2():分别实现一维信号和二维信号的 DWT（离散小波变换）idwt(),idwt2():分别实现一维信号和二维信号的 IDWT（离散小波变换） ；Wavedec2():多级二维小波分解函数；Waveinfo():提供小波包中所有的小波信息；攻击函数对算法进行攻击测试是对水印鲁棒性检测的一种重要手段，一个好的水印算法必须经过各种攻击测试才能对之做出客观的评价。MATLAB 中的许多函数可以直接用来做攻击测试。旋转：rotate（）可以对图像进行任意角度的旋
16、转；剪裁：imcrop()可以按精确定位的各点坐标进行剪裁；滤波：filter()和 filter2()可实现对一维信号和二维信号的滤波；抖动：dither()对图像进行抖动；抖动攻击考验水印鲁棒性的一个很好的攻击；jpeg 压缩：imwrite()中 jpg 和 quality 参数能对图像进行可控 jpg 压缩；加各种噪声：imnoise()可以对图像加入各种噪声，如白噪声、椒盐噪声等，加入噪声是对水印鲁棒性考验的一种常见的攻击；放大/缩小：imresize()可以以指定的插值方法来对图像进行放大和缩小。第四章： 图像数字水印技术的实现4.1 基于离散余弦变法（DCT）实现数字水印技术 打
17、开原始及水印图像：subplot(2,2,1)I=uigetfile(*.bmp,打开原始彩色图像文件 );RGB=imread(I);image(RGB);title(原始彩色图像);subplot(2,2,2)I=uigetfile(*.bmp,打开水印灰度图像文件 );imshow(I);title(灰度水印图像);subplot(2,2,3)H=imread(I);J=dct2(H);imshow(log(abs(J),),colorbar;title(水印图像经 DCT 变换后能量分布情况)运行结果： 水印全过程：0%水印加入程序Q=input(请输入放缩因子的值( 建议小于 1)：
18、Q )subplot(2,3,1)RGB=imread(浙江台州学院,jpg);imshow(RGB);title(原始图像);subplot(2,3,2)N=dct2(RGB(:,:,3);imshow(log(abs(N),),colorbar;title(Y 分量能量分布);subplot(2,3,4)I=imread(lena1,bmp);imshow(I);title(灰度水印图像);subplot(2,3,5)M=dct2(I);imshow(log(abs(M),),colorbar;title(水印能量分布);subplot(2,3,6)J=M(1:128,1:128);J(1
19、28:364,128:400)=0;J=rot90(J);J=rot90(J);J(365:600,401:750)=0;J=rot90(J);J=rot90(J);N=N+Q*J;K=idct2(N);RGB(:,:,3)=K;imshow(RGB);title(加入水印后图像);%水印提取程序subplot(2,3,3)RGB1=imread(浙江台州学院,jpg);N=dct2(RGB(:,:,3);M=dct2(RGB1(:,:,3);M=(N-M)/Q;B=idct2(M(236:365,350:401);Y=mat2gray(B);imshow(Y);title(提取的水印图像)运
20、行结果： 水印全过程（经剪切检测水印）%水印加入程序Q=input(请输入放缩因子的值( 建议小于 1)：Q )subplot(3,3,1)RGB=imread(MM,jpg);imshow(RGB);title(原始图像);subplot(3,3,2)imshow(RGB(:,:,3);title(B 分量);subplot(3,3,3)N=dct2(RGB(:,:,3);imshow(log(abs(N),),colorbar;title(B 分量能量分布);subplot(3,3,4)I=imread(lena1,bmp);imshow(I);title(灰度水印图像);subplot(
21、3,3,5)M=dct2(I);imshow(log(abs(M),),colorbar;title(水印能量分布);subplot(3,3,7)J=M(1:128,1:128);J(128:464,128:364)=0;J=rot90(J);J=rot90(J);J(465:800,365:600)=0;J=rot90(J);J=rot90(J);N=N+Q*J;K=idct2(N);RGB(:,:,3)=K;imshow(RGB);title(加入水印后图像);subplot(3,3,8)I=imcrop(RGB,1 1 598 798);imshow(I);subplot(3,3,9)%
22、水印提取程序subplot(3,3,6)RGB1=imread(MM,jpg);J=RGB1(:,:,3);X=J(1:799,1:599);N=dct2(I(:,:,3);M=dct2(X);M=(N-M)/Q;B=idct2(M(337:464,237:364);Y=mat2gray(B);imshow(Y);title(经放缩后提取的水印图像)运行结果： 水印全过程（经空域压缩检测水印）程序源代码%水印加入程序Q=input(请输入放缩因子的值( 建议小于 1)：Q )P=input(请输入您所希望的图像放缩系数值( 建议取值不要小于 0.5):P=)subplot(3,3,1)RGB=
23、imread(浙江台州学院,jpg);imshow(RGB);title(原始图像);subplot(3,3,2)imshow(RGB(:,:,3);title(B 分量);subplot(3,3,3)N=dct2(RGB(:,:,3);imshow(log(abs(N),),colorbar;title(B 分量能量分布);subplot(3,3,4)I=imread(lena1,bmp);imshow(I);title(灰度水印图像);subplot(3,3,5)M=dct2(I);imshow(log(abs(M),),colorbar;title(水印能量分布);subplot(3,3
24、,7)J=M(1:128,1:128);J(128:364,128:400)=0;J=rot90(J);J=rot90(J);J(365:600,401:750)=0;J=rot90(J);J=rot90(J);N=N+Q*J;K=idct2(N);RGB(:,:,3)=K;imshow(RGB);title(加入水印后图像);subplot(3,3,8)I=imresize(RGB,P,nearest);imshow(I);title(压缩 P 倍图像);subplot(3,3,9)J=imresize(I,1/P,nearest);imshow(J);title(再放大 P 倍还原图像)%
25、水印提取程序subplot(3,3,6)RGB1=imread(浙江台州学院,jpg);N=dct2(J(:,:,3);M=dct2(RGB1(:,:,3);M=(N-M)/Q;B=idct2(M(236:365,350:401);Y=mat2gray(B);imshow(Y);title(经放缩后提取的水印图像)：运行结果：4.2 图像水印的 dwt 算法%以下是水印提取算法clear all;clc;%保存时间start_time=cputime;figure(1);%读出原始图像subplot(1,2,1);input=imread(2.jpg);imshow(input);title(
26、原始图像);%读出水印图像subplot(1,2,2);watermarked_image=imread(watermarked.bmp);imshow(watermarked_image,);title(水印图像);%三色分离input=double(input);watermarked_image=double(watermarked_image);inputr=input(:,:,1);watermarked_imager=watermarked_image(:,:,1);inputg=input(:,:,2);watermarked_imageg=watermarked_image(:
27、,:,2);inputb=input(:,:,3);watermarked_imageb=watermarked_image(:,:,3);%水印图像R的分解Cwr,Swr=WAVEDEC2(watermarked_imager,2,haar);%图像R的分解Cr,Sr=WAVEDEC2(inputr,2,haar);%水印图像G的分解Cwg,Swg=WAVEDEC2(watermarked_imageg,2,haar);%图像R的分解Cg,Sg=WAVEDEC2(inputg,2,haar);%水印图像B的分解Cwb,Swb=WAVEDEC2(watermarked_imageb,2,haa
28、r);%图像B的分解Cb,Sb=WAVEDEC2(inputb,2,haar);%提取水印小波系数%提取水印R的小波系数r=0.06;for k=0:3whr(k+1,:)=Cwr(1+size(Cwr,2)/4+k*size(Cwr,2)/16:.size(Cwr,2)/4+(k+1)*size(Cwr,2)/16)-.Cr(1+size(Cr,2)/4+k*size(Cr,2)/16:.size(Cr,2)/4+(k+1)*size(Cr,2)/16);wvr(k+1,:)=Cwr(1+size(Cwr,2)/2+k*size(Cwr,2)/16:.size(Cwr,2)/2+(k+1)*
29、size(Cwr,2)/16)-.Cr(1+size(Cr,2)/2+k*size(Cr,2)/16:.size(Cr,2)/2+(k+1)*size(Cr,2)/16);wdr(k+1,:)=Cwr(1+3*size(Cwr,2)/4+k*size(Cwr,2)/16:.3*size(Cwr,2)/4+(k+1)*size(Cwr,2)/16)-.Cr(1+3*size(Cr,2)/4+k*size(Cr,2)/16:.3*size(Cr,2)/4+(k+1)*size(Cr,2)/16);endwhr=(whr(1,:)+whr(2,:)+whr(3,:)+whr(4,:)/(4*r);w
30、vr=(wvr(1,:)+wvr(2,:)+wvr(3,:)+wvr(4,:)/(4*r);wdr=(wdr(1,:)+wdr(2,:)+wdr(3,:)+wdr(4,:)/(4*r);war=(Cwr(1:size(Cwr,2)/16)-Cr(1:size(Cr,2)/16)/r;%提取水印小波系数%提取水印G的小波系数g=0.03;for k=0:3whg(k+1,:)=Cwg(1+size(Cwg,2)/4+k*size(Cwg,2)/16:.size(Cwg,2)/4+(k+1)*size(Cwg,2)/16)-.Cg(1+size(Cg,2)/4+k*size(Cg,2)/16:.s
31、ize(Cg,2)/4+(k+1)*size(Cg,2)/16);wvg(k+1,:)=Cwg(1+size(Cwg,2)/2+k*size(Cwg,2)/16:.size(Cwg,2)/2+(k+1)*size(Cwg,2)/16)-.Cg(1+size(Cg,2)/2+k*size(Cg,2)/16:.size(Cg,2)/2+(k+1)*size(Cg,2)/16);wdg(k+1,:)=Cwg(1+3*size(Cwg,2)/4+k*size(Cwg,2)/16:.3*size(Cwg,2)/4+(k+1)*size(Cwg,2)/16)-.Cg(1+3*size(Cg,2)/4+k*
32、size(Cg,2)/16:.3*size(Cg,2)/4+(k+1)*size(Cg,2)/16);endwhg=(whg(1,:)+whg(2,:)+whg(3,:)+whg(4,:)/(4*g);wvg=(wvg(1,:)+wvg(2,:)+wvg(3,:)+wvg(4,:)/(4*g);wdg=(wdg(1,:)+wdg(2,:)+wdg(3,:)+wdg(4)/(4*g);wag=(Cwg(1:size(Cwg,2)/16)-Cg(1:size(Cg,2)/16)/g;%提取水印小波系数%提取水印B的小波系数b=0.12;for k=0:3whb(k+1,:)=Cwb(1+size(
33、Cwb,2)/4+k*size(Cwb,2)/16:.size(Cwb,2)/4+(k+1)*size(Cwb,2)/16)-.Cb(1+size(Cb,2)/4+k*size(Cb,2)/16:.size(Cb,2)/4+(k+1)*size(Cb,2)/16);wvb(k+1,:)=Cwb(1+size(Cwb,2)/2+k*size(Cwb,2)/16:.size(Cwb,2)/2+(k+1)*size(Cwb,2)/16)-.Cb(1+size(Cb,2)/2+k*size(Cb,2)/16:.size(Cb,2)/2+(k+1)*size(Cb,2)/16);wdb(k+1,:)=C
34、wb(1+3*size(Cwb,2)/4+k*size(Cwb,2)/16:.3*size(Cwb,2)/4+(k+1)*size(Cwb,2)/16)-.Cb(1+3*size(Cb,2)/4+k*size(Cb,2)/16:.3*size(Cb,2)/4+(k+1)*size(Cb,2)/16);endwhb=(whb(1,:)+whb(2,:)+whb(3,:)+whb(4,:)/(4*b);wvb=(wvb(1,:)+wvb(2,:)+wvb(3,:)+wvb(4,:)/(4*b);wdb=(wdb(1,:)+wdb(2,:)+wdb(3,:)+wdb(4,:)/(4*b);wab=(
35、Cwb(1:size(Cwb,2)/16)-Cb(1:size(Cb,2)/16)/b;%重构水印图像cwr=war,whr,wvr,wdr;swr(:,1)=sqrt(size(war,2),sqrt(size(war,2),2*sqrt(size(war,2);swr(:,2)=sqrt(size(war,2),sqrt(size(war,2),2*sqrt(size(war,2);wr = waverec2(cwr,swr,haar);cwg=wag,whg,wvg,wdg;swg(:,1)=sqrt(size(wag,2),sqrt(size(wag,2),2*sqrt(size(wa
36、g,2);swg(:,2)=sqrt(size(wag,2),sqrt(size(wag,2),2*sqrt(size(wag,2);wg=waverec2(cwg,swg,haar);cwb=wab,whb,wvb,wdb;swb(:,1)=sqrt(size(wab,2),sqrt(size(wab,2),2*sqrt(size(wab,2);swb(:,2)=sqrt(size(wab,2),sqrt(size(wab,2),2*sqrt(size(wab,2);wb=waverec2(cwb,swb,haar);%将R，G，B叠加temp=size(wr);pic=zeros(temp
37、(1),temp(2),3);for i=1:temp(1);for j=1:temp(2);pic(i,j,1)=wr(i,j);pic(i,j,2)=wg(i,j);pic(i,j,3)=wb(i,j);endendoutput=uint8(round(pic);%转化为uint8watermark_image_uint8=uint8(output);imwrite(watermark_image_uint8,watermark.bmp,bmp);figure(2);imshow(watermark_image_uint8);title(提取出的水印);原始图像加入小波的图像小波 提取的小波
