一、实验目的

1. 了解数字图像水印技术的基本原理、分类和应用。
2. 掌握简单的可见水印和不可见水印的嵌入方法。
3. 实现一种基于DCT的不可见鲁棒水印，并进行水印鲁棒性测试。

二、实验内容

1. 简单的可见水印嵌入实验

图1给出了原始Lena图像和水印图像。请完成以下实验：
(1) 根据教材P445例8.29，将水印图像嵌入Lena图像(图像嵌入位置不唯一)，嵌入强度 $\alpha =3$。给出嵌入水印后的图像以及原始图像与嵌入水印图像的差值图像。(为更好地展示差值图像，在输出时可用函数imshow(filename, []))
(2) 请分别给出嵌入强度 $\alpha$ 为0、0.5 和 1 时的含水印图像，并分析 $\alpha$ 的改变对水印图像以及含水印图像主观质量的影响。

答
（1）

clear all;close all;clc
f = imread('1.Lena图像.tif');
wat_mark = imread('1.水印图像.tif');
a = 0.3;
pic_ori = f(200:251, 103:255); 
water_mark_pic = (1-a).*pic_ori + a.*wat_mark;
f_water_marked = f;
f_water_marked(200:251, 103:255) = water_mark_pic;
figure;imshow(f_water_marked);

（2）

dif=double(f)-double(f_water_marked);
figure;imshow(dif,[])%将差值归一化到（0~255）之间

%%

a=[0,0.5,1];

for i=1:3
    alpha = a(i);
    pic_ori = f(200:251, 103:255);
    water_mark_pic = (1-alpha).*pic_ori + alpha.*wat_mark;
    f_water_marked = f;
    f_water_marked(200:251, 103:255) = water_mark_pic;
    subplot(1,3,i);imshow(f_water_marked);
end

2. 不可见脆弱水印实验

(1) 请根据式(1)将水印图像(图1(b))嵌入到Lena图像(图1(a))的最低两位有效位(Least Significant Bit, LSB)中，并给出嵌入水印后的图像。(图像嵌入位置不唯一)
$$f_w=4\lfloor \frac{f}{4}\rfloor +\lfloor \frac{w}{64}\rfloor$$

其中，$f$ 为原始图像，$w$ 为水印图像，$f_w$ 为嵌入水印后的图像， 为下取整操作(可用函数floor)。

(2) 从含水印图像 $f_w$ 中提取水印图像 $w$。

(3) 对含水印图像进行 JPEG 压缩(质量因子Q=80)，以测试水印的脆弱性。给出压缩后的含水印图像以及提取的水印图像，并根据教材P389式(8.10)计算压缩后的含水印图像与原始含水印图像的均方根误差 $e_{\mathrm{rms}}$ (Root Mean Square Error, RMSE) 。

答

（1）

%在最低两位有效位嵌入水印

clear all;close all;clc
f = imread('1.Lena图像.tif');
m = imread('1.水印图像.tif');
f=double(f);%不使用double类型提取的水印效果较差
m=double(m);
W = f(204:255, 103:255); 
w=4.*floor(W./4)+floor(m./64);
fm = f;
fm(204:255, 103:255) = w;
fm=uint8(fm);
figure;imshow(fm)

（2）

%尝试提取水印
wat_mask=mod(fm,4);
wat_mask=wat_mask*64;
wat_mask=uint8(wat_mask);
figure;imshow(wat_mask,[]);

（3）

%对含水印图像进行JPEG压缩(质量因子Q=80)

imwrite(fm, '1.压缩图像.jpg', 'quality', 80)
c = imread('1.压缩图像.jpg'); 
C = double(c);


wat_mark_ext=c(204:255, 103:255);%水印提取

%次序从右向左排列，最右边的二进制位称为“最低有效位”
for i = 3:8
    wat_mark_ext = bitset(wat_mark_ext, i, 0); 
end

figure, imshow(wat_mark_ext,[])

%均方根
Fm = double(fm);%原始含水印图像
[r,c]=size(C);
erms=sqrt(sum(sum((Fm-C).^2))/(r*c))

erms =
2.3480

3. 不可见鲁棒水印实验

图2给出了原始Lena图像。请完成以下实验：
(1) 根据教材P447例8.30，实现基于DCT的不可见鲁棒水印。将一组伪随机序列作为水印，嵌入到图像较大的DCT系数中，获得含水印图像。嵌入强度 $\alpha =0.1$，伪随机序列长度 $k=1000$。
(2) 从含水印图像中提取水印，并计算原始水印与提取水印的相关系数 $\gamma$，利用二值检测判决(教材P449式(8.73))检测水印是否存在(阈值T设为0.8230)。
(3) 对含水印图像进行攻击测试，包括JPEG压缩(质量因子Q=80和20)、中值滤波(窗口大小3×3)、高斯噪声(均值0.1、方差0.001)、直方图均衡和图像旋转( $0.5^{\circ }$)。给出各种攻击下的含水印图像并计算原始水印与提取水印的相关系数 。

答

（1）

%将一组伪随机序列作为水印，嵌入到图像较大的DCT系数中，
% 获得含水印图像。嵌入强度0.1
% clear all;close all;clc
f = imread('3.Lena图像.tif');

f = double(f);
[r,c]=size(f);
fdct = dct2(f);
%对于固定的key，不是说从此以后产生的随机数都是相同的，
%而是在相同的key下，第一次调用rand产生的结果是相同的。
randn('seed', 0);
K = randn(1, 1000);
%将所有 DCT 系数变为一维行向量，然后按照幅值进行降序排序
fdct=reshape(fdct, 1, r*c);

[D,I] = sort(fdct,'descend');

%嵌入强度,选择长度为1000的伪随机序列作为无意义水印嵌入到图像中
a=0.1;
for i = 1:1000
    fdct(I(i)) = D(i)*(1+a*K(i));
end
%还原
fdct = reshape(fdct, r,c);

f_watermarked=idct2(fdct);
f_watermarked=uint8(f_watermarked);
figure,imshow(f_watermarked)



%从含水印图像中提取水印
% 并计算原始水印与提取水印的相关系数 gamma

a=0.1;T=0.8230;
w_hat=get_pseudo_random_sequence(f_watermarked,f,a,I);%得到估计的伪随机序列
[gamma,ifwatermarking]=water_extract(w_hat,T);%计算得到相关系数
disp(ifwatermarking);

（2）（3）

%%

%JPEG压缩攻击
imwrite(f_watermarked, '压缩图像.jpg', 'quality', 80);
f_JPEG_compressed =imread('压缩图像.jpg');
w_hat=get_pseudo_random_sequence(f_JPEG_compressed,f,a,I);%得到估计的伪随机序列
[gamma,ifwatermarking]=water_extract(w_hat,T);%计算得到相关系数
disp(ifwatermarking);

imwrite(f_watermarked, '压缩图像.jpg', 'quality', 20);
f_JPEG_compressed =imread('压缩图像.jpg');
w_hat=get_pseudo_random_sequence(f_JPEG_compressed,f,a,I);%得到估计的伪随机序列
[gamma,ifwatermarking]=water_extract(w_hat,T);%计算得到相关系数
disp(ifwatermarking);

%%
%平均滤波器攻击
% 定义平均核
kernel = ones(5, 5) / 25;
% 对每个通道应用平滑核
f_avg_filtered = imfilter(f_watermarked, kernel, 'replicate');
imshow(f_avg_filtered);
w_hat=get_pseudo_random_sequence(f_avg_filtered,f,a,I);%得到估计的伪随机序列
[gamma,ifwatermarking]=water_extract(w_hat,T);%计算得到相关系数
disp(ifwatermarking);

%%
%高斯噪声攻击，均值为0.1，方差为0.001的噪声
[M,N]=size(f);
gaussian_noise=imnoise2('gaussian',M,N,0.1,0.001);
f_gau_noise=f+gaussian_noise;
w_hat=get_pseudo_random_sequence(f_gau_noise,f,a,I);%得到估计的伪随机序列
[gamma,ifwatermarking]=water_extract(w_hat,T);%计算得到相关系数
disp(ifwatermarking);
%%
%中值滤波器攻击
f_median_filtered= medfilt2(f_watermarked, [3, 3]);
imshow(f_median_filtered);
w_hat=get_pseudo_random_sequence(f_median_filtered,f,a,I);%得到估计的伪随机序列
[gamma,ifwatermarking]=water_extract(w_hat,T);%计算得到相关系数
disp(ifwatermarking);

%%
%直方图均衡攻击
f_histogram_equalization= imadjust(f_watermarked);
imshow(f_histogram_equalization);
w_hat=get_pseudo_random_sequence(f_histogram_equalization,f,a,I);%得到估计的伪随机序列
[gamma,ifwatermarking]=water_extract(w_hat,T);%计算得到相关系数
disp(ifwatermarking);
%%
%图像旋转0.5度攻击
imshow(f_watermarked);
f_rotated=imrotate(f_watermarked,0.5,'bicubic', 'crop');
figure;imshow(f_rotated);
w_hat=get_pseudo_random_sequence(f_rotated,f,a,I);%得到估计的伪随机序列
[gamma,ifwatermarking]=water_extract(w_hat,T);%计算得到相关系数
disp(ifwatermarking);
%%

function w_hat=get_pseudo_random_sequence(f_watermarked,f,a,I)
f_watermarked=double(f_watermarked);
f_watermarked_dct=dct2(f_watermarked);

f_watermarked_dct=reshape(f_watermarked_dct, 1, []);

%根据初始的图像f生成原始的fdct
f = double(f);
[r,c]=size(f);
fdct = dct2(f);
fdct=reshape(fdct, 1, r*c);

w_hat=zeros(1,1000);
for i = 1:1000
    w_hat(i) = (f_watermarked_dct(I(i)) - fdct(I(i))) / (a * fdct(I(i)));
end
end


function [gamma,ifwatermarking] = water_extract(w_hat,T)
randn('seed', 0);
w = randn(1, 1000);
w_avg=mean(w);
w_hat=double(w_hat);
w_hat_avg=mean(w_hat);

X1 = w - w_avg;
Y1 = w_hat - w_hat_avg;
num = sum(X1.*Y1);
deno = sqrt(sum(X1.^2)*sum(Y1.^2));
gamma=num/deno;

if gamma>=T
    ifwatermarking='图像含有水印';
else
    ifwatermarking='图像不含有水印';
end
end