实验2：基于LSB和DCT的图像数字水印

1. 实验类别

        设计型实验：MATLAB设计并实现基于LSB的图像数字水印算法、基于DCT的图像数字 水印算法，计算峰值信噪比值。

2. 实验目的

        了解信息隐藏中最常用的LSB算法的特点，掌握LSB算法原理，设计并实现一种基于图像的LSB隐藏算法。

        了解基于DCT的图像数字水印技术，掌握基于DCT系数关系的图像水印算法原理，设计并实现一种基于DCT的数字水印算法。

        了解如何通过峰值信噪比来对图像数字水印效果进行客观评价。

3. 实验条件

  • Windows 2000或Windows Xp以上操作系统；
  • MATLAB 6.5以上版本软件；
  • 图像文件

4. 实验原理

点击查看

4.1 基于LSB的图像数字水印

        任何多媒体信息在数字化时都会产生物理随机噪声，而人的感官系统对这些随机噪声并不敏感。替换技术就是利用这个原理，通过使用秘密信息比特替换随机噪声，从而实现信息隐藏目的。图像高位平面对图像感官质量起主要作用，去除图像最低几个位平面并不会造成画面质量的下降。利用这个原理可用秘密信息（或称水印信息）替代载体图像低位平面以实现信息嵌入。
        LSB算法选用最低位平面来嵌入秘密信息，最低位平面对图像的视觉效果影响最轻微，但很容易受噪声影响和攻击，可采用冗余嵌入的方式来增强稳健性加以解决，即在一个区域中嵌入相同的信息，提取时根据该区域中的所有像素判断。

4.2 基于DCT的图像数字水印

        在信号的频域（变换域）中隐藏信息要比在时域中嵌入信息具有更好的鲁棒性。一幅图像经过时域到频域的变换后，可将待隐藏信息嵌入图像的显著区域，这种算法更具抗攻击能力，而且保持了人类感官的不可察觉性。常用的变换域方法有离散余弦变换、离散小波变换和离散傅立叶变换等。
        本实验介绍一种提取秘密信息的时候不需要原始图像的盲水印算法，算法的思想是利用载体中两个特定的DCT系数的相对大小来表示隐藏的信息。首先载体图像分为8*8分块，进行二维DCT变换，分别选择其中约定的两个位置，比如用（u1,v1）和（u2,v2）代表所选定的两个系数的坐标。如果Bi（u1,v1）< Bi（u2,v2），代表隐藏1。如果Bi（u1,v1）> Bi（u2,v2），代表隐藏0。
        提取的时候接收者对包含水印的图像文件进行二维DCT变换，比较每一块中约定位置的DCT系数值，根据其相对大小，如果Bi（u1,v1）< Bi（u2,v2），则提取1；如果Bi（u1,v1）> Bi（u2,v2），则提取0。最终得到隐藏信息的比特串，从而恢复出秘密信息。
        在上述算法过程中，如果有一对系数大小相差非常少，往往难以保证携带图像在保存和传输的过程中以及提取秘密信息的过程中不发生变化。因此在实际的设计过程中，一般都是引入一个Alpha变量对系数的差值进行控制，将两个系数的差值放大，可以保证提取秘密信息的正确性。

4.3 峰值信噪比

        峰值信噪比（常缩写为PSNR）是一个表示信号最大可能功率和影响它的表示精度的破坏性噪声功率的比值的工程术语。由于许多信号都有非常宽的动态范围，峰值信噪比常用对数分贝单位来表示。峰值信噪比经常用作图像压缩等领域中信号重建质量的测量方法，它常简单地通过均方误差（MSE）进行定义。

5. 实验步骤

5.1.基于LSB的图像数字水印

第一步：隐藏算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
clear;clc;
msgfid = fopen('hidden.txt','r'); % 打开秘密文件，读入秘密信息
[msg,count] = fread(msgfid,'ubit1'); % 读取秘密信息存入msg，count为成功读入了多少位
fclose(msgfid); % 关闭文件
img = imread('lena.bmp'); % 读入图像
img_R = img(:,:,1); % 提取RGB第1层嵌入水印
[row,col] = size(img_R); % 获取img_R有row行col列
max_count = row * col; % 图像能嵌入水印最大比特数

% 判断是否将秘密信息可以完全嵌入进去
if count > max_count
    disp('图片容量无法通过LSB方法嵌入所有水印信息!');quit;
else
    disp(['图片通过LSB能嵌入的最大水印比特数为: ',num2str(max_count),' bits']);
end

msg_counter = 1; % 记录秘密信息存入进度
total_watermark_bits = 0; % 记录图片写入多少位秘密信息bit

% 计算需要在图片中嵌入多少行，mod取余，floor向下取整
if mod(count,row) == 0
    round = floor(count/row);
else
    round = floor(count/row) + 1;
end

% 遍历嵌入秘密信息，bitset设置指定位置的位
for outer=1:1:round
    for inner = 1:1:col
       img_R(inner,outer) = bitset(img_R(inner,outer),1,msg(msg_counter,1));% 改变像素值最后一位bit为水印bit值
       msg_counter = msg_counter + 1;
       total_watermark_bits = total_watermark_bits + 1;
       if msg_counter > count % 嵌入完毕跳出内循环
           break; 
       end
    end  % 内层
    if msg_counter > count % 嵌入完毕跳出外循环
        disp(['LSB嵌入正常结束!共嵌入水印比特数: ',num2str(total_watermark_bits),' bits']);
        break;
    end
end  % 外层

img(:,:,1) = img_R; % 覆盖
imwrite(img,'./LSB_watermarked.bmp'); % 保存为水印图像
figure('NumberTitle', 'off', 'Name', '图像水印：LSB');
subplot(1,2,1);imshow('lena.bmp');title('原始图像');
subplot(1,2,2);imshow(img);title('嵌入水印图像');

运行结果如图，可以看出原始图像与嵌入水印图像看不出明显差别。

第二步：提取算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
clear;clc;
watermark_bits = 202960; % 水印比特数
fid = fopen('LSB_message.txt','w');   %打开输出秘密的文件
img = imread('LSB_watermarked.bmp'); % 读取含有水印信息的彩色图像
img_R = img(:,:,1); % 提取RGB第1层嵌入水印
msg = zeros(watermark_bits,1); % 创建一个watermark_bits行1列的全0矩阵，用于存放水印比特
[row,col] = size(img_R);% 获取img_R有row行col列

msg_counter = 1; % 记录秘密信息提取进度
total_watermark_bits = 0; % 记录图片提取了多少位信息bit

% 计算图片中的水印信息存了多少行，mod取余，floor向下取整
if mod(watermark_bits,row) == 0
    round = floor(watermark_bits/row);
else
    round = floor(watermark_bits/row) + 1;
end

% 遍历提取秘密信息，bitset设置指定位置的位
for outer = 1:1:round
    for inner = 1:1:col
        msg(msg_counter,1) = bitget(img_R(inner,outer),1);  % 提取图像矩阵的bit水印信息
        msg_counter = msg_counter + 1;
        total_watermark_bits = total_watermark_bits + 1;
        if msg_counter > watermark_bits % 提取完毕跳出内循环
            break;
        end 
    end % 内层
    if msg_counter > watermark_bits
        disp(['LSB水印提取正常结束!共提取水印比特数:',num2str(total_watermark_bits),' bits']);
        break; % 提取完毕跳出内循环
    end
end % 外层

fwrite(fid,msg,'ubit1'); % 存入提取的信息，以ubit1的格式
fclose(fid); % 关闭文件

5.2.基于DCT的图像数字水印

第一步：隐藏算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
clear;clc;
msgfid = fopen('hidden.txt','r'); % 打开秘密文件，读入秘密信息
[msg,count] = fread(msgfid,'ubit1'); %读取秘密信息存入msg，count为成功读入了多少位
fclose(msgfid); % 关闭文件

I = imread('lena.bmp'); % 读入图像
I = double(I)/255; % 将图片转化为double型
I_R= I(:,:,1); % 取图象的R分量
T = @dct2; 
DCTrgb = blkproc(I_R,[8 8],T); % 图像分块变换
[row,col] = size(DCTrgb); % 获取图片的row行col列
row = floor(row/8); % 分块后的行数
col = floor(col/8); % 分块后的列数

% 判断是否超出容量
if(count > row * col) 
    disp('水印内容超过图片容量');quit;
end

k = 1; % 记录嵌入秘密信息比特数
temp = 0; % 交换的中间变量
alpha=0.02; % 设置alpha变量

% 选择(5,2)和(4,3)这一对系数，对分块后的这两个位置进行调整从而实现嵌入水印
for i = 0:row-1
    for j = 0:col-1
        if msg(k,1) == 0
            if DCTrgb(i*8+5,j*8+2) < DCTrgb(i*8+4,j*8+3) 
                temp = DCTrgb(i*8+5,j*8+2);
                DCTrgb(i*8+5,j*8+2) = DCTrgb(i*8+4,j*8+3);
                DCTrgb(i*8+4,j*8+3) = temp;
            end 
        else
            if DCTrgb(i*8+5,j*8+2) > DCTrgb(i*8+4,j*8+3)
                temp = DCTrgb(i*8+5,j*8+2);
                DCTrgb(i*8+5,j*8+2) = DCTrgb(i*8+4,j*8+3);
                DCTrgb(i*8+4,j*8+3) = temp;
            end
        end
        % 放大两个系数的差值，把原本小的系数调整到更小，即把（5,2）和（4,3）中的最小值再减去alpha；
        if DCTrgb(i*8+5,j*8+2) < DCTrgb(i*8+4,j*8+3)
            DCTrgb(i*8+5,j*8+2) = DCTrgb(i*8+5,j*8+2) - alpha;
        else 
            DCTrgb(i*8+4,j*8+3) = DCTrgb(i*8+4,j*8+3) - alpha;
        end
        k = k+1;
        if k > count
            disp(['嵌入秘密文件信息的比特数: ',num2str(k-1),' bits']);
            break; % 内层
        end
    end
    if k > count
        break; %外层
    end
end
T = @idct2;
wi = blkproc(DCTrgb,[8 8],T);      %将DCT进行逆变换
I(:,:,1) = wi; % 把更改后的R图层替换回原来彩色图像的R层
imwrite(I,'DCT_watermarked.bmp');  %保存为水印图像
figure('NumberTitle', 'off', 'Name', '图像水印：DCT');
subplot(1,2,1);imshow('lena.bmp');title('原始图像');
subplot(1,2,2);imshow('DCT_watermarked.bmp');title('嵌入水印图像');

第二步：提取算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
clear;clc;
count = 704; % 秘密信息的比特数
fid = fopen('DCT_message.txt','w'); % 打开输出秘密的文件
wi = imread('DCT_watermarked.bmp'); % 读取含有水印信息的彩色图像
wi = double(wi)/255; % 将图片转化为double型
wi_R = wi(:,:,1); % 取图象的R分量
T = @dct2; 
DCTcheck = blkproc(wi_R,[8 8],T); % 对图像分块进行DCT变换
[row,col] = size(DCTcheck); % 获取图片的row行col列
row = floor(row/8); % 分块后的行数
col = floor(col/8); % 分块后的列数

k = 1; % 记录已提取秘密信息的比特数
% 判断(5,2),(4,3)的系数值，其实就是嵌入的逆过程
for i = 0:row-1
    for j = 0:col-1
        if DCTcheck(i*8+5,j*8+2) < DCTcheck(i*8+4,j*8+3)
            fwrite(fid,1,'ubit1'); % 写入1
        else
            fwrite(fid,0,'ubit1'); % 写入0
        end
        if k > count
            break; % 内层
        end
        k = k+1;
    end
    if k > count
        break; % 外层
    end
end
fclose(fid);  % 关闭文件

5.3.峰值信噪比计算

$MAX = 2^8-1$

$MSE(均方误差) = \frac{1}{mn}\sum{i=0}^{m-1}\sum{j=0}^{n-1}[I(i,j)-K(i,j)]^2$

$PSNR(峰值信噪比)=10·log{10}(\frac{MAX_I^2}{MSE}) = 20·log{10}(\frac{MAX_I}{\sqrt{MSE}})$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
clear;clc;
W = imread('lena.bmp');
W = W(:,:,1);

E = imread('LSB_watermarked.bmp'); % 以LSB水印为例
E = E(:,:,1);

F = W - E;
imshow(double(F));
[m,n] = size(W);

MAX = 2^8 - 1;
MSE = sum(sum((W-E).^2))/(m*n);
apsnr = 20*log10(MAX/sqrt(MSE));
disp('--------LSB峰值信噪比计算--------')
disp(['MAX = ',num2str(MAX)]);
disp(['MSE = ',num2str(MSE)]);
disp(['PSNR = ',num2str(apsnr)]);

6. 实验报告

1．理解整个数字水印嵌入和提取过程，试写出实验的关键步骤。简单概述你在编程过程中遇到了哪些问题，如何解决的？

  • 实验内容及步骤如上所示。
  • 学习认识blkproc函数，其为图像进行分块处理函数，通过参数可以在分块的同时进行图像变换。
  • blkproc和blockproc：MATLAB分块处理矩阵的函数很早就有了——blkproc，但缺点是blkproc必须一次性把要处理的矩阵全部导入内存中，这样大大限制了其应用范围，对一些超大型的图像就无能为力了。随着MATLAB使用范围越来越广，Mathworks也紧跟用户需求，新版本的MATLAB推出了可以处理任意大图像的函数blockproc。但是在使用时，将blkproc替换为blockproc后出现报错，老老实实用blkproc吧。
  • matlab中@的使用:用于定义函数句柄的操作符。函数句柄既是一种变量，可以用于传参和赋值；也是可以当做函数名一样使用。
  • DCT解密出后发现是乱码，排查出示嵌入水印代码的问题，debug半天发现是其索引i和j的错误，无语。。。

2．通过对嵌入水印后图像的攻击，分析该算法的鲁棒性和不可察觉性如何，Alpha变量对鲁棒性和不可察觉性又有什么影响？

        a.通过噪声，滤波，压缩，剪切，旋转分别对LSB和DCT算法下的水印进行攻击，分析其峰值信噪比来判断其鲁棒性和不可察觉性。（数字水印算法设计有五种攻击的实验过程及结论，此处只以旋转攻击为例）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
% DCT旋转攻击，LSB同理
watermarked = imread('DCT_watermarked.bmp'); % 读取DCT算法水印图像
W = watermarked; E = watermarked;
E = imrotate(E,10,'bilinear','crop'); % 对水印图像进行旋转10度
imwrite(E,'DCT旋转攻击.bmp') % 存储图像，用于DCT提取。

% 峰值信噪比计算
F = W - E;
imshow(double(F));
[m,n] = size(W);

MAX = 2^8 - 1;
MSE = sum(sum((W-E).^2))/(m*n);
apsnr = 20*log10(MAX/sqrt(MSE));
disp('--------DCT(a=0.02)旋转攻击峰值信噪比计算--------')
disp(['MAX = ',num2str(MAX)]);
disp(['MSE = ',num2str(MSE)]);
disp(['PSNR = ',num2str(apsnr)]);
figure('NumberTitle', 'off', 'Name', 'DCT(a=0.02)旋转攻击');
subplot(1,2,1);imshow(W);title('原始图像');
subplot(1,2,2);imshow(E);title('旋转10度');

        如上图所示的LSB和DCT在旋转攻击下，PNSR明显降低，旋转后的图片提取水印为乱码，说明两个水印算法抗旋转攻击能力弱。根据参考文献可知，在抗噪声，滤波，压缩攻击有很好的能力，但是对于抗裁剪和旋转攻击能力较弱。

        b.更改DCT算法下的Alpha值分别为0.01，0.02，0.05，0.08时的PNSR，从而研究Alpha变量对鲁棒性和不可察觉性的影响。

        在水印嵌入过程中，为了满足不可见性的条件，就要要求水印嵌入的强度不能太高，但是为了能在受攻击后准确地恢复出水印信息，通常期望嵌入的水印强度越大越好，因此要在不可见性和鲁棒性之间寻找一个折中。根据上图实验所示，随着alpha的增大，PSNR值减小，不可见性变差。由参考文献图片所示，alpha越大，相似度系数NC越接近原图像。因此alpha要取中频，以做到鲁棒性和不可见性都较满足。

7.可选实验

        顺序LSB算法是顺序选取像素点进行信息隐秘，会造成图像各部分统计特征的不一致，而导致了严重的安全问题，修改部分与未修改部分会具有不同的统计特性，增大了攻击者对秘密通信怀疑的可能性。为了解决这个问题，可以采用随机间隔选取像素序列，进行随机修改嵌入，试编写基于密钥控制的随机lsb隐写算法。

        因为水印内容过多，会导致随机数过长不好输入测试，因此水印内容取：L。过程中发现字符的比特是小端存储。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
% 随机间隔LSB嵌入
clear;clc;
msgfid = fopen('hidden.txt','r'); % 打开秘密文件，读入秘密信息
[msg,count] = fread(msgfid,'ubit1'); % 读取秘密信息存入msg，count为成功读入了多少位
fclose(msgfid); % 关闭文件
img = imread('lena.bmp'); % 读入图像
img_R = img(:,:,1); % 提取RGB第1层嵌入水印
[row,col] = size(img_R); % 获取img_R有row行col列
max_count = row * col; % 图像能嵌入水印最大比特数

% 产生随机序列，即为密钥
flag_1 = 0; % 记录随机序列有多少个1，满足水印的比特数
flag_size = 0;
while flag_1 < count
    disp('-----正在产生随机序列-----');
    % 产生可以覆盖图像的随机0,1序列
    random = randi([0 1],1,max_count);
    % 循环判断其中的1是否满足水印所需的比特数
    for i = 1:max_count
        flag_size = flag_size+1; % 要放在循环前面，break后会少算
        if random(i) == 1
            flag_1 = flag_1 + 1;
            if flag_1 >= count
                break;
            end
        end
    end
end
random = random(1:flag_size); % 保留有效位，删除冗余的随机数
disp('随机序列产生完毕：');
disp(random);
    
msg_counter = 1; % 记录秘密信息存入进度
total_watermark_bits = 0; % 记录图片写入多少位秘密信息bit

% 计算需要在图片中嵌入多少行，mod取余，floor向下取整
if mod(flag_size,row) == 0
    round = floor(flag_size/row);
else
    round = floor(flag_size/row) + 1;
end

% 遍历嵌入秘密信息，bitset设置指定位置的位
for outer=1:1:round
    for inner = 1:1:col
       total_watermark_bits = total_watermark_bits + 1;
       if random(1,outer * col + inner -col) == 1 
           img_R(inner,outer) = bitset(img_R(inner,outer),1,msg(msg_counter,1));% 改变像素值最后一位bit为水印bit值
           msg_counter = msg_counter + 1;
           if msg_counter > count % 嵌入完毕跳出内循环
               break; 
           end
       end
    end % 内层
    if msg_counter > count % 嵌入完毕跳出外循环
        disp(['LSB嵌入正常结束!共嵌入水印比特数: ',num2str(total_watermark_bits),' bits']);
        break;
    end
end  % 外层

img(:,:,1) = img_R; % 覆盖
imwrite(img,'./LSB_watermarked.bmp'); % 保存为水印图像
figure('NumberTitle', 'off', 'Name', '图像水印：LSB随机嵌入');
subplot(1,2,1);imshow('lena.bmp');title('原始图像');
subplot(1,2,2);imshow(img);title('嵌入随机水印图像');

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
clear;clc;
random_key = [1,0,1,1,1,0,1,0,0,0,0,1,0,1,1]; % 嵌入时的随机序列
watermark_bits = 14; % 水印比特数
fid = fopen('LSB_message.txt','w');   %打开输出秘密的文件
img = imread('LSB_watermarked.bmp'); % 读取含有水印信息的彩色图像
img_R = img(:,:,1); % 提取RGB第1层嵌入水印
msg = zeros(watermark_bits,1); % 创建一个watermark_bits行1列的全0矩阵，用于存放水印比特
[row,col] = size(img_R);% 获取img_R有row行col列

msg_counter = 1; % 记录秘密信息提取进度
total_watermark_bits = 0; % 记录图片提取了多少位信息bit

% 计算图片中的水印信息存了多少行，mod取余，floor向下取整
if mod(watermark_bits,row) == 0
    round = floor(watermark_bits/row);
else
    round = floor(watermark_bits/row) + 1;
end

% 遍历提取秘密信息，bitset设置指定位置的位
for outer = 1:1:round
    for inner = 1:1:col
        total_watermark_bits = total_watermark_bits + 1;
        if random_key(1,outer * col  - col + inner) == 1
            msg(msg_counter,1) = bitget(img_R(inner,outer),1);  % 提取图像矩阵的bit水印信息
            msg_counter = msg_counter + 1;
            if total_watermark_bits >= watermark_bits % 提取完毕跳出内循环,需要等于不然密钥数组越界
                break;
            end
        end
    end % 内层
    if total_watermark_bits >= watermark_bits
        disp(['LSB水印提取正常结束!共提取水印比特数:',num2str(total_watermark_bits),' bits']);
   ,break; % 提取完毕跳出内循环
    end
end % 外层

msg = msg(1:msg_counter-1); % 删除一开始0矩阵没用用到的空间
fwrite(fid,msg,'ubit1'); % 存入提取的信息，以ubit1的格式
fclose(fid); % 关闭文件

鸣谢❀参考大佬文章

  • 基于LSB的图像数字水印实验
  • 数字水印算法设计
  • LSB顺序+随机隐写和提取(matlab)
  • 陈铭. 隐写与隐写分析算法及实践研究[D].北京邮电大学,2008.
  • 林洪文,杨绍清. 基于DCT变换的数字水印算法[A]. 中国通信学会.第十八届全国青年通信学术年会论文集（上册）[C].中国通信学会:中国通信学会,2013:4.
  • 李倩,宋弘. 基于图像数字处理在DCT域的应用[A]. 中国人工智能学会智能数字内容安全专业委员会.Proceedings of 2010 Asia-Pacific Conference on Information Network and Digital Content Security（2010APCID）[C].中国人工智能学会智能数字内容安全专业委员会:中国通信学会青年工作委员会,2010:6.