我在研究如何对图像进行二维快速傅里叶变换，求指点

zola · 公众号 · 自媒体 · 2017-02-23 23:34

正文

我发现了一个给图片加不可见的水印的方法，据说这种“盲水印”不影响用户体验，又能追溯图片来源和传播过程，在捍卫版权和知识付费越来越看重的现在，这个东西应该蛮有前途的。我就在想，我没有Matlab2010a这软件，我无法体验给图片加密，我试试找找能否用PHP实现在线加盲水印的办法，希望能做成WEB应用。

我得到的线索来自这里：

解密：阿里巴巴公司根据截图查到泄露信息的员工的技术是什么

https://www.zhihu.com/question/50735753

解密：阿里巴巴公司根据截图查到泄露信息的员工的技术是什么

http://www.phpchina.com/portal.php?mod=view&aid=40242

大概技术细节是这两张图来解释：

加水印：

然后解水印：

里面关键的部分是对图片进行“二维快速傅里叶变换”，英文叫 Two-Dimensional Fast Fourier Transform，在技术文档里通常简称为FFT。

我阅读了Fuqiang Liu 的Matlab2010a下的源代码，也去读了《如果看了此文你还不懂傅里叶变换，那就过来掐死我吧【完整版】》，初步了解了快速傅里叶变换是怎么回事。

我找到以下关于FFT的PHP代码实现相关链接：

https://www.phpclasses.org/package/6193-PHP-Compute-the-Fast-Fourier-Transform-of-sampled-data.html ，这个类库能Compute the Fast Fourier Transform of sampled data，但我不知道怎么使用，有样本代码教我怎么做FFT和inverse FFT。

http://www.jasonbailey.net/stuff/php-fast-fourier-transform-fft-brighton-php-october-2013/

我也尝试找PHP是否支持快速傅里叶变换，找到有图像组件

http://php.net/manual/zh/book.imagick.php 下的 Imagick::forwardFourierTransformImage 是支持 “离散傅里叶变换”，英文是discrete Fourier transform (DFT) ，但不知道与“快速傅里叶变换”Fast Fourier Transform（FFT）有什么差别。

这里还有一个 Two-Dimensional Fast Fourier Transform（二维快速傅里叶变换）的文章，

http://imagejdocu.tudor.lu/doku.php?id=gui:process:fft 也是基于 Matlab在讲。

现在问题是，用php里的imagecopy函数合并背景图片和水印到同一张照片的功能我测试通过了，如何用FFT的类库对图片作FFT和inverse FFT操作？求有图像处理经验的的高手指导一下，我已经试过用英文搜索了，如 php Fourier image / php fft convert image / ,找到这些链接，http://www.fmwconcepts.com/misc_tests/FFT_tests/

http://www.roborealm.com/help/FFT.php

都不知道如何使用傅里叶变换对图片进行操作。

下面是Fuqiang Liu 的Matlab2010a下的源代码：

%%傅里叶变换加水印源代码%% 运行环境Matlab2010a clc;clear;close all;alpha = 1;%% read dataim = double(imread('gl1.jpg'))/255;mark = double(imread('watermark.jpg'))/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');%imwrite(uint8(GG),'extracted watermark.jpg');%% MSE and PSNRC=double(im);RC=double(FAO);MSE=0; PSNR=0;for i=1:imsize(1)
    for j=1:imsize(2)
        MSE=MSE+(C(i,j)-RC(i,j)).^2;
    endendMSE=MSE/360.^2;PSNR=20*log10(255/sqrt(MSE));MSEPSNR%% attack test%% attack by smearing%A = double(imread('gl1.jpg'));%B = double(imread('attacked image.jpg'));attack = 1-double(imread('attack.jpg'))/255;figure,imshow(attack);FAO_ = FAO;for i=1:imsize(1)
    for j=1:imsize(2)
        if attack(i,j,1)+attack(i,j,2)+attack(i,j,3)>0.5
            FAO_(i,j,:) = attack(i,j,:);
        end
    endendfigure,imshow(FAO_);%extract watermarkFA2=fft2(FAO_);G=(FA2-FA)*2;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');%% attack by cuttings2 = 0.8;FAO_ = FAO;FAO_(:,s2*imsize(2)+1:imsize(2),:) = FAO_(:,1:int32((1-s2)*imsize(2)),:);figure,imshow(FAO_);%extract watermarkFA2=fft2(FAO_);G=(FA2-FA)*2;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');%%小波变换加水印，解水印大家按照加的思路逆过来就好clc;clear;close all;%% read dataim = double(imread('gl1.jpg'))/255;mark = double(imread('watermark.jpg'))/255;figure, imshow(im),title('original image');figure, imshow(mark),title('watermark');%% RGB divisionim=double(im); mark=double(mark); imr=im(:,:,1); markr=mark(:,:,1); img=im(:,:,2); markg=mark(:,:,2); imb=im(:,:,3); markb=mark(:,:,3); %% parameterr=0.04; g = 0.04; b = 0.04;%% wavelet tranform and add watermark% for red[Cwr,Swr]=wavedec2(markr,1,'haar'); [Cr,Sr]=wavedec2(imr,2,'haar'); % add watermarkCr(1:size(Cwr,2)/16)=... Cr(1:size(Cwr,2)/16)+r*Cwr(1:size(Cwr,2)/16); k=0; while k=size(Cr,2)/size(Cwr,2)-1 Cr(1+size(Cr,2)/4+k*size(Cwr,2)/4:size(Cr,2)/4+... (k+1)*size(Cwr,2)/4)=Cr(1+size(Cr,2)/4+... k*size(Cwr,2)/4:size(Cr,2)/4+(k+1)*size(Cwr,2)/4)+... r*Cwr(1+size(Cwr,2)/4:size(Cwr,2)/2); Cr(1+size(Cr,2)/2+k*size(Cwr,2)/4:size(Cr,2)/2+... (k+1)*size(Cwr,2)/4)=Cr(1+size(Cr,2)/2+... k*size(Cwr,2)/4:size(Cr,2)/2+(k+1)*size(Cwr,2)/4)+... r*Cwr(1+size(Cwr,2)/2:3*size(Cwr,2)/4); Cr(1+3*size(Cwr,2)/4+k*size(Cwr,2)/4:3*size(Cwr,2)/4+... (k+1)*size(Cwr,2)/4)=Cr(1+3*size(Cr,2)/4+... k*size(Cwr,2)/4:3*size(Cr,2)/4+(k+1)*size(Cwr,2)/4)+... r*Cwr(1+3*size(Cwr,2)/4:size(Cwr,2)); k=k+1; end; Cr(1:size(Cwr,2)/4)=Cr(1:size(Cwr,2)/4)+r*Cwr(1:size(Cwr,2)/4); % for green[Cwg,Swg]=WAVEDEC2(markg,1,'haar'); [Cg,Sg]=WAVEDEC2(img,2,'haar'); Cg(1:size(Cwg,2)/16)=... Cg(1:size(Cwg,2)/16)+g*Cwg(1:size(Cwg,2)/16); k=0; while k=size(Cg,2)/size(Cwg,2)-1 Cg(1+size(Cg,2)/4+k*size(Cwg,2)/4:size(Cg,2)/4+... (k+1)*size(Cwg,2)/4)=Cg(1+size(Cg,2)/4+... k*size(Cwg,2)/4:size(Cg,2)/4+(k+1)*size(Cwg,2)/4)+... g*Cwg(1+size(Cwg,2)/4:size(Cwg,2)/2); Cg(1+size(Cg,2)/2+k*size(Cwg,2)/4:size(Cg,2)/2+... (k+1)*size(Cwg,2)/4)=Cg(1+size(Cg,2)/2+... k*size(Cwg,2)/4:size(Cg,2)/2+(k+1)*size(Cwg,2)/4)+... g*Cwg(1+size(Cwg,2)/2:3*size(Cwg,2)/4); Cg(1+3*size(Cg,2)/4+k*size(Cwg,2)/4:3*size(Cg,2)/4+... (k+1)*size(Cwg,2)/4)=Cg(1+3*size(Cg,2)/4+... k*size(Cwg,2)/4:3*size(Cg,2)/4+(k+1)*size(Cwg,2)/4)+... g*Cwg(1+3*size(Cwg,2)/4:size(Cwg,2)); k=k+1; end; Cg(1:size(Cwg,2)/4)=Cg(1:size(Cwg,2)/4)+g*Cwg(1:size(Cwg,2)/4); % for blue[Cwb,Swb]=WAVEDEC2(markb,1,'haar'); [Cb,Sb]=WAVEDEC2(imb,2,'haar'); Cb(1:size(Cwb,2)/16)+b*Cwb(1:size(Cwb,2)/16); k=0; while k=size(Cb,2)/size(Cwb,2)-1 Cb(1+size(Cb,2)/4+k*size(Cwb,2)/4:size(Cb,2)/4+... (k+1)*size(Cwb,2)/4)=Cb(1+size(Cb,2)/4+... k*size(Cwb,2)/4:size(Cb,2)/4+(k+1)*size(Cwb,2)/4)+... g*Cwb(1+size(Cwb,2)/4:size(Cwb,2)/2); Cb(1+size(Cb,2)/2+k*size(Cwb,2)/4:size(Cb,2)/2+... (k+1)*size(Cwb,2)/4)=Cb(1+size(Cb,2)/2+... k*size(Cwb,2)/4:size(Cb,2)/2+(k+1)*size(Cwb,2)/4)+... b*Cwb(1+size(Cwb,2)/2:3*size(Cwb,2)/4); Cb(1+3*size(Cb,2)/4+k*size(Cwb,2)/4:3*size(Cb,2)/4+... (k+1)*size(Cwb,2)/4)=Cb(1+3*size(Cb,2)/4+... k*size(Cwb,2)/4:3*size(Cb,2)/4+(k+1)*size(Cwb,2)/4)+... b*Cwb(1+3*size(Cwb,2)/4:size(Cwb,2)); k=k+1; end; Cb(1:size(Cwb,2)/4)=Cb(1:size(Cwb,2)/4)+b*Cwb(1:size(Cwb,2)/4); %% image reconstructionimr=WAVEREC2(Cr,Sr,'haar'); img=WAVEREC2(Cg,Sg,'haar'); imb=WAVEREC2(Cb,Sb,'haar'); imsize=size(imr); FAO=zeros(imsize(1),imsize(2),3); for i=1:imsize(1); for j=1:imsize(2); FAO(i,j,1)=imr(i,j); FAO(i,j,2)=img(i,j); FAO(i,j,3)=imb(i,j); end end figure, imshow(FAO); title('watermarked image');

地