任务背景及简介

 随着互联网的普及和数字技术的广泛应用，信息媒体的数字化极大地提高了信息的存取效率, Intenet 的普及又为数字信息在世界范围内的迅速传递开辟了便捷的途径，数字产品变得极大丰富并且传播便捷，但与此同时版权保护问题也日益突出。传统的信息安全技术在数字产品版权保护方面存在的不足，促成了数字水印技术的发展。数字水印技术将数字水印隐藏于数字媒体中,以便在版权纠纷中为版权所有者提供版权证明；然而数字水印技术处在起步阶段,许多技术还不成熟，很多问题和方法有待进一步的研宄和解决。

本实验课题主要是研究数字水印在数字图像中的嵌入、检测提取的方法，在水印算法中引入图像的变换域用于嵌入水印，并利用人类的视觉系统对不同颜色的敏感程度不同以此增强水印的透明性。

算法原理及实现

理论基础

（1）DCT变换

在水印算法中引入图像的变换域用于嵌入水印，以此增强水印的透明性。（将一组光强数据转换成频率数据，以便得知强度变化的情形。若对高频的数据做些修饰，再转回原来形式的数据时，显然与原始数据有些差异，但是人类的眼睛却是不容易辨认出来。）

（2）人类视觉系统特性

人类的视觉系统对不同颜色的敏感程度不同。人眼对绿色光最为敏感，对红色光的敏感程度次之，对于蓝色光最不敏感。（人眼对红色光和蓝色光的敏感程度之和与对绿色光的敏感程度较为接近。）该实验课题使用红色和蓝色光部分作为 Pacthwork 算法的 A 集合，绿色光作为 B 集合，两个集合在嵌入信息时使用相逆的操作，可以在一定程度上互相抵消嵌入信息引起的图像视觉上的变化，提高水印透明性。

具体算法

（1）水印的嵌入

图片预处理：

载体图像original与水印图像mark选取正方形的RGB图像，并满足载体图像的边长为 8 的倍数且载体图像边长为水印图像边长的8倍。将载体图像的 3个颜色通道分离，得到IＲ、IG、IB 3个颜色分量; 将水印图像的3个颜色通道分离，得到WＲ、WG、WB 3个颜色分量。

对水印Arnold置换：

对WR、WG、WB三个颜色分量进行Arnold变换，（变换的过程可以看作是拉伸、压缩、折叠及拼接的过程）经过变换得到WRA、WGA、WBA。Arnold变换过程如下：

其中(x,y)为图像上点坐标。

对载体图像进行DCT变换并提取直流分量：

载体图像的各分量以8×8的大小为一个单位划分为若干子块。将子块视为一个整体，最左上角的子块位置坐标为(1,1)，其相邻的右边的子块位置坐标为(1，2) ，依此类推。对每一个子块分别应用DCT 变换，然后取出变换后的每一个子块左上角的直流分量组成一个新矩阵，位置坐标为(1,1) 的子块的直流分量作为新矩阵( 1,1) 位置的元素，位置坐标为 (1,2) 的子块的直流分量作为新矩阵(1,2) 位置的元素，依此类推。最终所得的矩阵称为直流分量矩阵IRD、IBD、IGD。

在直流分量矩阵中嵌入水印：

在直流分量矩阵上嵌入水印，嵌入方法是增加/减去置乱的水印图像分量k倍的亮度。分量提取公式如下：

其中绿色分量在亮度处理时使用减操作，红色分量和蓝色分量使用加操作。（由人类视觉系统的特性知这两个互逆的操作可以在一定程度上互相抵消嵌入信息引起的图像视觉上的变化，提高水印透明性。）

合成嵌入水印图像：

嵌入水印后的直流分量矩阵 IRDE、IBDE、IGDE按照对应位置替换各个子块的直流分量，再对各个子块分别应用反DCT 变换完成各颜色分量的水印嵌入。

最后将三个分量合成为水印图像。

（2）水印的提取

预处理：

将含有数字水印的载体图像I_embed的3个颜色通道分离，得到PR、PG、PB 3个颜色分量; 将原始载体图像original的3个颜色通道分离，得到IR、IG、IB 3 个颜色分量。

对含有水印的图像进行分块 DCT 变换并提取直流分量：

对带水印载体图像P使用和水印嵌入时相同直流分量提取方式，得到PRD、PBD、PGD。

从直流分量矩阵中提取置乱的水印：

将三个分量进行逆Arnold变换后组合成提取出的水印。

实验结果及分析

嵌入水印并提取：

已知原图original为jpg格式，水印mark为png格式，直接对图片进行加水印和提取水印操作

直观上来看，可以发现嵌入载体后和嵌入载体前图像肉眼无法分辨差异，但是提取出的水印图像和原水印间存在一定失真。

结果分析：

（1）隐蔽性分析：

利用python脚本对载体图像和水印图像进行互相关系数计算：

从计算结果可以看出两张图片具有99%的相似性，该算法具有隐蔽性。

（2）容量分析：从容量上来看，假设该算法载体图像大小为mm，水印图像可为(m/8)(m/8)

（3）鲁棒性检验：对嵌入水印的图片I_embed进行白噪声攻击、高斯滤波攻击、裁剪攻击和旋转攻击

白噪声攻击

对进行攻击后的水印与原水印进行互相关系数计算

2、高斯滤波攻击

计算互相关系数

3、裁剪攻击：

计算互相关系数

4、旋转攻击：

计算互相关系数

从上述各实验效果可看出该算法具有一定鲁棒性。

参考博客：https://blog.csdn.net/qq_44465615/article/details/116420445?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522162320718416780262565873%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fall.%2522%257D&request_id=162320718416780262565873&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~first_rank_v2~rank_v29-7-116420445.pc_search_result_control_group&utm_term=matlab%E5%AE%9E%E7%8E%B0%E6%95%B0%E5%AD%97%E6%B0%B4%E5%8D%B0%E7%AE%97%E6%B3%95%EF%BC%8C%E9%AA%8C%E8%AF%81%E9%B2%81%E6%A3%92%E6%80%A7&spm=1018.2226.3001.4187

计算图片的互相关系数（python脚本）
from PIL import Image
from numpy import average, linalg, dot
 
 
def get_thumbnail(image, size=(1200, 750), greyscale=False):
    image = image.resize(size, Image.ANTIALIAS)
    if greyscale:
        image = image.convert('L')
    return image
def image_similarity_vectors_via_numpy(image1, image2):
    image1 = get_thumbnail(image1)
    image2 = get_thumbnail(image2)
    images = [image1, image2]
    vectors = []
    norms = []
    for image in images:
        vector = []
        for pixel_tuple in image.getdata():
            vector.append(average(pixel_tuple))
        vectors.append(vector)
        norms.append(linalg.norm(vector, 2))
    a, b = vectors
    a_norm, b_norm = norms
    res = dot(a / a_norm, b / b_norm)
    return res
image1 = Image.open('F:\TEST4\mark.png')
image2 = Image.open('F:\TEST4\whirl_mark.png')
cosin = image_similarity_vectors_via_numpy(image1, image2)
print(cosin)

%% 初始化
clc;
clear ;
figure(1);      %打开窗口
%% 载入数据
I=imread("original.jpg");
%I=imresize(I,[512,512],'nearest'); 
subplot( 2,2,1) ,imshow(I),title('载体图像');
W=imread("mark.png");
%W=imresize(W,[64,64],'nearest');
subplot(2,2,2),imshow(W),title('水印图像');
[row,col,t] = size(W);
 
%% 嵌入水印部分
%分离R，G，B通道
IR=I(:,:,1);
IG=I(:,:,2);
IB=I(:,:,3);
WR=W(:,:,1);
WG=W(:,:,2);
WB=W(:,:,3);
%设置k,不同图片k值不同
k=0.162;
 
%% 使用Arnold变换置乱水印
WRA=arnold(WR,1,1,1);
WRA=double(WRA);
WGA=arnold(WG,1,1,1);
WGA=double(WGA);
WBA=arnold(WB,1,1,1);
WBA=double(WBA);
%%
% 对载体图像进行8*8分块处理，然后对每块分别DCT变化
IRD=blkproc(IR,[8,8],'dct2');
IGD=blkproc(IG,[8,8],'dct2');
IBD=blkproc(IB,[8,8],'dct2');
IRDE=IRD;
IGDE=IGD;
IBDE=IBD;
 
%% 提取直流分量,并向矩阵中嵌入水印
for i=0:(row-1)
    for j=0:(col-1)
        x=i*8;
        y=j*8;
        IRDE(x+1,y+1)=IRD(x+1,y+1)+k*WRA(i+1,j+1);
        IBDE(x+1,y+1)=IBD(x+1,y+1)+k*WBA(i+1,j+1);
        IGDE(x+1,y+1)=IGD(x+1,y+1)-k*WGA(i+1,j+1);
    end
end
 
%% 对载体图像进行分块反DCT变换
IR2=blkproc(IRDE,[8,8],'idct2');
IG2=blkproc(IGDE,[8,8],'idct2');
IB2=blkproc(IBDE,[8,8],'idct2');
IR2=uint8(IR2);
IG2=uint8(IG2);
IB2=uint8(IB2);
 
 
%% 合成
I_embed = I;
I_embed(:,:,1) = IR2;
I_embed(:,:,2) = IG2;
I_embed(:,:,3) = IB2;
subplot( 223) ,imshow(I_embed),title('嵌入水印后的载体图像');
 
%% 提取水印
%% 分离通道
P=I_embed;
PR=P(:,:,1);
PG=P(:,:,2);
PB=P(:,:,3);
 
%% 将带水印图像进行DCT变换
PRD=blkproc(PR,[8,8],'dct2');
PGD=blkproc(PG,[8,8],'dct2');
PBD=blkproc(PB,[8,8],'dct2');
WR2=WR;
WB2=WB;
WG2=WG;
%% 提取水印
for i=0:(row-1)
    for j=0:(col-1)
        x=i*8;
        y=j*8;
        WR2(i+1,j+1)=(PRD(x+1,y+1)-IRD(x+1,y+1))/k;
        WB2(i+1,j+1)=(PBD(x+1,y+1)-IBD(x+1,y+1))/k;
        WG2(i+1,j+1)=(IGD(x+1,y+1)-PGD(x+1,y+1))/k;
    end
end
%% 逆arnold 
WR2=uint8(WR2);
WG2=uint8(WG2);
WB2=uint8(WB2);
WR2=rearnold(WR2,1,1,1);
WG2=rearnold(WG2,1,1,1);
WB2=rearnold(WB2,1,1,1);
%% 合成水印
W2=W;
W2(:,:,1)=WR2;
W2(:,:,2)=WG2;
W2(:,:,3)=WB2;
subplot( 224) ,imshow(W2),title('提取出的水印');
imwrite(I_embed,'I_embed.png');
imwrite(W2,'W2.png');
%% arnold变换
function arnoldImg = arnold(img,a,b,n)
[h,w] = size(img);
N=h;
arnoldImg = zeros(h,w);
for i=1:n
    for y=1:h
        for x=1:w
            %防止取余过程中出现错误，先把坐标系变换成从0 到 N-1
            xx=mod((x-1)+b*(y-1),N)+1;
            yy=mod(a*(x-1)+(a*b+1)*(y-1),N)+1;  
            arnoldImg(yy,xx)=img(y,x);              
        end
    end
    img=arnoldImg;
end
arnoldImg = uint8(arnoldImg);
end
 
function img = rearnold(arnoldImg,a,b,n)
[h,w] = size(arnoldImg);
img = zeros(h,w);
N = h;
for i=1:n
    for y=1:h
        for x=1:w           
            xx=mod((a*b+1)*(x-1)-b*(y-1),N)+1;
            yy=mod(-a*(x-1)+(y-1),N)+1  ;      
            img(yy,xx)=arnoldImg(y,x);              
        end
    end
    arnoldImg=img;
end