【新星计划】Matlab-傅里叶变换隐藏水印

看到一个文章，关于阿里巴巴公司根据截图查到泄露信息的具体员工的技术是什么，这个新闻有点老了，主要说一下背后的图片隐藏水印方法

一、水印原理

相对于空域方法，频域加盲水印的方法隐匿性更强，抵抗攻击能力更强。这类算法解水印困难，你不知道水印加在那个频段，而且受到攻击往往会破坏图像原本内容。

所谓盲水印，是指人感知不到的水印，包括看不到或听不见（没错，数字盲水印也能够用于音频）。其主要应用于音像作品、数字图书等，目的是，在不破坏原始作品的情况下，实现版权的防护与追踪。

添加数字盲水印的方法简单可分为空域方法和频域方法，这两种方法添加了冗余信息，但在编码和压缩情况不变的情况下，不会使原始图像大小产生变化（原来是10MB添加盲水印之后还是10MB）。

空域是指空间域，我们日常所见的图像就是空域。空域添加数字水印的方法是在空间域直接对图像操作（之所以说的这么绕，是因为不仅仅原图是空域，原图的差分等等也是空域），比如将水印（Matlab）直接叠加在图像上。

我们常说一个音有多高，这个音高是指频率；同样，图像灰度变化强烈的情况，也可以视为图像的频率。频域添加数字水印的方法，是指通过某种变换手段（傅里叶变换，离散余弦变换，小波变换等）将图像变换到频域（小波域），在频域对图像添加水印，再通过逆变换，将图像转换为空间域。相对于空域手段，频域手段隐匿性更强，抗攻击性更高。

所谓对水印的攻击，是指破坏水印，包括涂抹，剪切，放缩，旋转，压缩，加噪，滤波等。数字盲水印不仅仅要敏捷性高（不被人抓到），也要防御性强（抗打）。就像Dota的敏捷英雄往往是脆皮，数字盲水印的隐匿性和鲁棒性是互斥的。（鲁棒性是抗攻击性的学术名字）

 

二、频域制作数字盲水印的方法
信号是有频率的，一个信号可以看做是无数个不同阶的正弦信号的的叠加。

上式为傅里叶变换公式，f(t)是指时域信号（对于信号我们说时域，因为是与时间有关的，而图像我们往往说空域，与空间有关），F(w)是指频率。想要对傅里叶变换有深入了解的同学，建议看一下《信号与系统》或者《数字信号处理》的教材，里面系统介绍了傅里叶变换、快速傅里叶变换、拉普拉斯变换、z变换等。

简而言之，我们有方法将时域信号转换成为频域，同样，我们也能将二维信号（图像）转换为频域。在上文中提到，图像的频率是指图像灰度变换的强烈情况。

下面以傅里叶变换为例，介绍通过频域给图像添加数字盲水印的方法。注意，因为图像是离散信号，我们实际用的是离散傅里叶变换，在本文采用的都是二维快速傅里叶变换，快速傅里叶变换与离散时间傅里叶变换等价，下文中傅里叶变换均为二维快速傅里叶变换。

 

上图为叠加数字盲水印的基本流程。编码的目的有二，一是对水印加密，二控制水印能量的分布。

 

三、叠加数字盲水印的实验

给出一张200*400的原图

之后进行傅里叶变换，下图变换后的频域图像，

水印图为60*200

编码后的水印，编码方式采用随机序列编码，通过编码，水印分布到随机分布到各个频率，并且对水印进行了加密，水印注意是对称的。

将上图与原图的频谱叠加，可见图像的频谱已经发生了巨大的变化

之后，将叠加水印的频谱进行傅里叶逆变换，得到叠加数字水印后的图像，

肉眼几乎看不出叠加水印后的图像与原图的差异，这样，数字盲水印已经叠加到图像中去。
实际上，我们是把水印以噪声的形式添加到原图像中。
下图是在空域上的加水印图与原图的残差

可以看出，实际上上述方法是通过频域添加冗余信息（像噪声一样）。这些噪声遍布全图，在空域上并不容易破坏。

那么，为什么频谱发生了巨大的变化，而在空域却变化如此小呢？这是因为我们避开了图像的主要频率。下图是原图频谱，其能量主要集中在低频。

水印提取是水印叠加的逆过程，

经提取后，我们得到如下水印，问：为什么水印要对称呢？

四、源代码

•  

clc;clear;close all;alpha = 1;
%% read dataim = double(imread('3.png'))/255;mark = double(imread('4.png'))/255;figure, imshow(im),title('original image');figure, imshow(mark),title('watermark');
%% encode markimsize = size(im);%randomTH=zeros(imsize(1)*0.5,imsize(2),imsize(3));TH1 = TH;TH1(1:size(mark,1),1:size(mark,2),:) = mark;M=randperm(0.5*imsize(1));N=randperm(imsize(2));save('encode.mat','M','N');for i=1:imsize(1)*0.5    for j=1:imsize(2)        TH(i,j,:)=TH1(M(i),N(j),:);    endend% symmetricmark_ = zeros(imsize(1),imsize(2),imsize(3));mark_(1:imsize(1)*0.5,1:imsize(2),:)=TH;for i=1:imsize(1)*0.5    for j=1:imsize(2)        mark_(imsize(1)+1-i,imsize(2)+1-j,:)=TH(i,j,:);    endendfigure,imshow(mark_),title('encoded watermark');%imwrite(mark_,'encoded watermark.jpg');
%% add watermarkFA=fft2(im);figure,imshow(FA);title('spectrum of original image');FB=FA+alpha*double(mark_);figure,imshow(FB); title('spectrum of watermarked image');FAO=ifft2(FB);figure,imshow(FAO); title('watermarked image');%imwrite(uint8(FAO),'watermarked image.jpg');RI = FAO-double(im);figure,imshow(uint8(RI)); title('residual');%imwrite(uint8(RI),'residual.jpg');xl = 1:imsize(2);yl = 1:imsize(1);[xx,yy] = meshgrid(xl,yl);figure, plot3(xx,yy,FA(:,:,1).^2+FA(:,:,2).^2+FA(:,:,3).^2),title('spectrum of original image');figure, plot3(xx,yy,FB(:,:,1).^2+FB(:,:,2).^2+FB(:,:,3).^2),title('spectrum of watermarked image');figure, plot3(xx,yy,FB(:,:,1).^2+FB(:,:,2).^2+FB(:,:,3).^2-FA(:,:,1).^2+FA(:,:,2).^2+FA(:,:,3).^2),title('spectrum of watermark');
%% extract watermarkFA2=fft2(FAO);G=(FA2-FA)/alpha;GG=G;for i=1:imsize(1)*0.5    for j=1:imsize(2)        GG(M(i),N(j),:)=G(i,j,:);    endendfor i=1:imsize(1)*0.5    for j=1:imsize(2)        GG(imsize(1)+1-i,imsize(2)+1-j,:)=GG(i,j,:);    endendfigure,imshow(GG);title('extracted watermark');