基于量化的图像水印算法

张道德 UCAS第2题 基于量化的图像水印算法抖动调制是一种特殊的QIM 方法，但结构简单，其特点是给定量化器的任何量化单元和重建点都是其他量化器量化单元和重建点的平移。

在无水印的情况下 ，平移量是随机的，成为抖动量。

为了嵌入信息，抖动量可以待嵌入的水印信息进行调制。

将各个可能的带嵌入的的水印信息映射成不同的抖动量()d W 。

载体信号在进行抖动后，通过量化产生合成信号。

在采用基本量化器()q ⋅的情况下，嵌入函数可以表示为(;) (())()E X W q round X d W d W =+-这类水印嵌入技术就称为抖动调制。

2.1 算法介绍下面介绍采取均匀标量量化器的二进制抖动调制算法：假设要嵌入的水印信息是一个二维序列，即{0,1}W =，其长度为L 。

根据DM 的嵌入函数，可以得到二值水印嵌入公式为()(()[,])[,]k k S k q X k d k b d k b =+-其中，()X k 是载体用来嵌入信息的的k 个数据，{0,1}k b ∈与水印W 相对应，[,]k d k b 是抖动量，()q ⋅是基本量化器，即()(/)q x round x ∆=∆⨯∆这里，∆是量化步长，x 是待量化数据，()round ⋅表示四舍五入操作，()q ∆⋅表示量化步长为∆的量化器。

• 确定要修改的载体数据()X k ，即从载体数据中根据一定原则选择合适的载体数据来嵌入水印信息；• 确定用于对各个载体数据()X k 进行量化的步长k ∆。

如果各个载体数据的量化步长不同，则是一种自适应量化；一般情况下，可以采用等步长的非适应量化，此时k ∆相等，均为∆；• 确定抖动量[,]k d k b ，一般[,0]d k 可任意选择，但是[,1]d k 要按下式进行选取：[,1][,0]/2, [,1]0[,1][,0]/2, [,1]0k k d k d k if d k d k d k if d k =+∆<⎧⎨=-∆>⎩在量化步长时，可以为每个载体数据选择相同的抖动量；• 根据公式()(()[,])[,]k k S k q X k d k b d k b =+-，先用量化器函数()q ⋅对载体数据进行量化，然后合成输出值()S k 。

在提取水印，一般采取最小距离检测，即{0,1}()arg min |[]|k k Y l b k Y S l ∧∈=- 其中，k Y 是检测器接收到的信号数据，[0]k Y S 和[1]k Y S 分别为假设k Y 是以[,0]d k 和[,1]d k为抖动量的情况下应该得到的抖动调制结果。

2.2编程实现与分析我们在本实验嵌入的水印信息为一张大小为3232⨯的二值黑白图像。

原始图像（载体图像）为“elain1.bmp ”，其大小为512512⨯。

我们将图像分块，分成88⨯的图像小块，每块图像中嵌入1比特水印数据信息，易知可以讲水印信息完整嵌入。

2.2.1抖动调制的实验结果与分析水印嵌入步骤如下：• 确定要修改的载体数据。

我们在本实验中选取两块图像小块的像素差值作为要修改的载体数据，位置为[5][5]和[4][4]，并将得到的信息保存在[4][4]，即[4][4][5][5](([5][5][4][4]))q d d ∆=--+-其中d 为抖动量。

• 我们取量化步长∆是固定的，即20∆=，是一种非适应量化。

• 确定抖动量[,]k d k b 。

我们在本实验中随机选择[,0][30,30]d k ∈-的整数， [,1]d k 按下式进行选取：[,1][,0]/2, [,1]0[,1][,0]/2, [,1]0k k d k d k if d k d k d k if d k =+∆<⎧⎨=-∆>⎩并将[,0]d k 保存在一个一位数组中，传给水印提取者。

按照上述方法，将水印信息（二值黑白图像）嵌入到载体图像当中，如下图所示：图2.1 原始图像与嵌入水印后的图像两张图像用肉眼观察没有区别，可以得出水印信息的嵌入没有对原始图像在视觉上造成水印提取采取最小距离检测从而提取水印。

计算如下数值：0011([5][5][4][4])(([5][5][4][4]))([5][5][4][4])(([5][5][4][4]))x q d dy q d d ∆∆=---+-⎧⎨=---+-⎩ 其中，0d 和1d 分别为当前嵌入信息为0和1的抖动量。

从而得到当前嵌入的水印信息0, 1, if x y b if x y∧<⎧=⎨<⎩ 提取水印结果很理想，与原始水印图像是一致的，结果如下图：图2.2 原始水印图像与提取的水印图像 将水印图像中值滤波后，再次提取水印图像。

对水印图像进行中值滤波后，与原始图像有些区别，中值滤波后的图像更加清晰，如下图：图2.3 原始图像与中值滤波后的水印图像 从中值滤波后的水印图像提取水印图像是很差的，几乎看不出原始水印图像的样子，如下图：图2.4 原始水印图像与中值滤波后提取的水印图像2.2.2 程序代码%************************************************************************************ ①求解水印图像的峰值信噪比和均方误差函数function[MseValue,PsnrValue]=MsePsnr(OrigImg,Watermark)%计算均方误差%OrigImg为原始图像%Watermark为嵌入水印信息后的图像MseValue=0;OrigImgMatrix=double(imread(OrigImg)); %读入更改后bmp位图,存入矩阵ImgMatrix中WatermarkMatrix=double(imread(Watermark));%读入更改后bmp位图,存入矩阵ImgMatrix中[r c]=size(OrigImgMatrix);for i=1:rfor j=1:cMseValue=MseValue+(OrigImgMatrix(i,j)-WatermarkMatrix(i,j))^2;endendMseValue=MseValue/(r*c);PsnrValue=10*log10(255^2/MseValue);%************************************************************************************* ②水印嵌入函数，抖动嵌入函数function[ArrayB,WMsize,ChangedMatrix]=QIM2(OrigImg,Watermark,step)WaterMatrix=double(imread(Watermark)); %读入水印图像,存入矩阵WaterMatrix中BiMtrix=im2bw(WaterMatrix); % 将水印图像转化为二值图像[Br,Bc]=size(BiMtrix);ImgMatrix=double(imread(OrigImg)); %读入bmp位图,存入矩阵ImgMatrix中[r,c]=size(ImgMatrix);ChangedMatrix=ImgMatrix; %更改后的图像数据r8=floor(r/8);c8=floor(c/8);BlockCell=cell(r8,c8);ArrayB=[]; %存储抖动量的数组WMsize=[Br,Bc]; %存储水印图像的大小% ********************************************************************************%将图像进行分块，每块为8*8的子块for i=1:r8for j=1:c8BlockCell{i,j}=ImgMatrix(((i-1)*8+1):(8*(i-1)+8),(8*(j-1)+1):(8*(j-1)+8));endend%对分块后的图像进行DCT变换for i=1:r8for j=1:c8BlockCell{i,j}=dct2(BlockCell{i,j});endend% ********************************************************************************%将水印图像嵌入到原始图像中for i=1:Brfor j=1:Bcd0=randi([-30 30],1); %抖动量d0ArrayB=[ArrayB d0];if d0<0d1=d0+step/2;elsed1=d0-step/2;endif BiMtrix(i,j)==0BlockCell{i,j}(4,4)=BlockCell{i,j}(5,5)-(step*round((BlockCell{i,j}(5,5)-BlockCell{i,j}(4,4)+d0)/ step)-d0);elseBlockCell{i,j}(4,4)=BlockCell{i,j}(5,5)-(step*round((BlockCell{i,j}(5,5)-BlockCell{i,j}(4,4)+d1)/ step)-d1);endendend% ********************************************************************************%执行逆ＤＣＴ逆变换，并将更改后的图像数据存储到矩阵ChangedMatrix中for i=1:r8for j=1:c8BlockCell{i,j}=idct2(BlockCell{i,j});ChangedMatrix(((i-1)*8+1):(8*(i-1)+8),(8*(j-1)+1):(8*(j-1)+8))=BlockCell{i,j};endend% ********************************************************************************③水印提取函数function [WaterMatrix]=ExtractWater(WMsize,ArrayB,ChangedImg,step)% 抖动调制算法，QIM的提取水印算法ImgMatrix=double(imread(ChangedImg)); %读入更改后bmp位图,存入矩阵ImgMatrix中[r,c]=size(ImgMatrix);r8=floor(r/8);c8=floor(c/8);BlockCell=cell(r8,c8);Br=WMsize(1); % 水印图像的行数Bc=WMsize(2);% 水印图像的列数WaterMatrix=zeros(WMsize);% ********************************************************************************%将图像进行分块，每块为8*8的子块for i=1:r8for j=1:c8BlockCell{i,j}=ImgMatrix(((i-1)*8+1):(8*(i-1)+8),(8*(j-1)+1):(8*(j-1)+8));endend%对分块后的图像进行DCT变换for i=1:r8for j=1:c8BlockCell{i,j}=dct2(BlockCell{i,j});endend% ********************************************************************************%将水印图像从更改后的图像中提取出来for i=1:Brfor j=1:Bcd0=ArrayB((i-1)*Bc+j); %抖动量if d0<0d1=d0+step/2;elsed1=d0-step/2;endx=(BlockCell{i,j}(5,5)-BlockCell{i,j}(4,4))-(step*round((BlockCell{i,j}(5,5)-BlockCell{i,j}(4,4)+d 0)/step)-d0);y=(BlockCell{i,j}(5,5)-BlockCell{i,j}(4,4))-(step*round((BlockCell{i,j}(5,5)-BlockCell{i,j}(4,4)+d 1)/step)-d1);if abs(x)<abs(y)WaterMatrix(i,j)=0;elseWaterMatrix(i,j)=1;endendend% ******************************************************************************** ④主函数clcclearstep=20; %量化步长%*********************************************************************** [ArrayB,WMsize,ChangedMatrix]=QIM2('elain1.bmp','watermark.pgm',step);imwrite(uint8(ChangedMatrix),'elain1_mark.bmp');[WaterMatrix]=ExtractWater(WMsize,ArrayB,'elain1_mark.bmp',step);imwrite(WaterMatrix,'watermark_Changed.pgm');figure(1)subplot(121),imshow('elain1.bmp'),title('原始图像');subplot(122),imshow('elain1_mark.bmp'),title('用QIM嵌入水印后的图像');figure(2)subplot(121),imshow('watermark.pgm'),title('原始水印图像');subplot(122),imshow('watermark_changed.pgm'),title('提取的水印图像');%***********************************************************************%对水印图像采取中值滤波imwrite(uint8(medfilt2(ChangedMatrix)),'elain1_mark2.bmp');[WaterMatrix2]=ExtractWater(WMsize,ArrayB,'elain1_mark2.bmp',step);imwrite(WaterMatrix2,'watermark_Changed2.pgm');figure(3)subplot(121),imshow('elain1.bmp'),title('原始图像');subplot(122),imshow('elain1_mark2.bmp'),title('中值滤波后的水印图像');figure(4)subplot(121),imshow('watermark.pgm'),title('原始水印图像');subplot(122),imshow('watermark_changed2.pgm'),title('中值滤波后提取的水印图像'); %***********************************************************************%计算水印图像的峰值信噪比与均方误差[MseValue,PsnrValue]=MsePsnr('elain1.bmp','elain1_mark.bmp');disp('水印图像的峰值信噪比与均方误差分别为')PsnrValueMseValue% ********************************************************************************。