【图像隐写】基于matlab LDPC编码译码改进DCT水印嵌入提取【含Matlab源码 832期】

一、DCT数字水印嵌入与提取简介

1 基本DCT变换
目前，基于DCT域的水印方法已经成为数字水印算法研究的热点，它的核心思想就是通过离散傅立叶变换对图像块进行处理后，再选择变换域中的一些系数值依据一定规则来嵌入水印。
由于图像块中DCT系数频带分布由左上角的直流分量DC往下对应的系数频率由低频升至高频，因此在不影响原图质量的前提下，可将水印信息根据能量大小嵌入相应系数频带中。通过图像块量化与水印嵌入结合的处理方法将水印信息均匀分布在图像的整个空间域，在图像裁剪和滤波方面，变换域的水印比在空间域的更能表现出一定的鲁棒性。

2 水印算法描述
2.1 水印嵌入算法
该算法采用加性嵌入的方式在经过DCT变换后的子图像块的中频域中，选取隐秘位置嵌入水印信息，具体的嵌入流程如下图1所示：

图1 分块水印嵌入流程
(1）分块处理：设宿主图像为P，将其分块处理为8*8的K个子块。
(2）水印预处理：设水印图像为W，对其进行互补变换，变换后的水印图像和变换前的水印图像相互补。
(3）对水印图像进行Arnold置乱变换，并依据混沌映射规则，选取密钥混沌序列并与水印序列异或运算，将置换次数和异或运算处理后的结果分别作为水印嵌入算法的密钥1和密钥2。
(4)DCT变换：对各子块内做DCT变换，利用zig-zag对DCT系数进行扫描，得到第k块子图像块的序列为Zk(i),i=0,1,2，…63.
(5）水印嵌入算法：依据zig-zag排序，在各子块的中、低频段选取特定系数x(m）和x(n），在系数坐标（a,b）和（c,d）处嵌入水印信息图像W，并将其作为密钥3。同理，嵌入互补水印图片W’，并将嵌入的位置作为密钥4。水印嵌入的方法如下：


(6)IDCT变换：将每一个子图像块作二维DCT逆变换。
(7）子块合并：将每一个子块合并成嵌入水印的图像P’。

2.2 水印提取算法
将嵌入水印的图像P’分块处理，并对各子块进行二维DCT变换，由密钥3和4推断所选择的水印系数，若x(m）≤x(n），则水印信息为0，若x(m)>x(n），则水印信息为1，再利用密钥1和2将初步水印的信息解密再进行Arnold逆变换，最终提取出水印信息。

2.3 水印检测算法
本文通过计算峰值信噪比PSNR的值评价嵌入水印的宿主图像的质量，一幅m和n的图像，PSNR度量标准定义为：

归一化相关系数NC的值判断嵌入水印的图像与宿主图像的相似度，其定义为：

二、部分源代码

close all;
clc;
clear all
warning off
%-----------------读入"隐藏的图片"---------------------
I=imread('W.bmp');
%-----------------------读入"载体图像"-------------------------
cover_image=imread('lena.bmp');
%------------------------------------------------------------------
I0=rgb2gray(I);%灰度化

cover_image=rgb2gray(cover_image);%灰度化
[wm0,watermarked_image,wm]=ldpc_dct(I0,cover_image);%ldpc_dct嵌入提取
e=wm0-wm;
[m,n]=size(e);
mse=sum((e(:).^2))/(m*n);
psnr=10*log10(255^2/mse);%原始水印与提取水印的峰值信噪比
% disp(['ldpc改进dct提取水印的峰值信噪比psnr=',num2str(psnr)])
figure(1)
subplot(221)
imshow(cover_image);
title('原图');
subplot(222);
imshow(I0);
title('水印图');
title('水印图');
%显示嵌入水印后的图象
subplot(223);
uint8_watermarked_image=uint8(watermarked_image);
imshow(uint8_watermarked_image)
title('ldpc编码译码改进后嵌入水印图')  
subplot(224);
imshow(double(wm));
title('ldpc编码译码改进后提取水印图')
%% 剪切攻击
I_jianqie=I0;%剪切图
I_jianqie(20:30,20:40)=256;
[wm_jianqie0,watermarked_image_jianqie,wm_jianqie]=ldpc_dct(I_jianqie,cover_image);%ldpc_dct嵌入提取
e_jianqie=wm_jianqie0-wm_jianqie;
[m,n]=size(e_jianqie);
mse_jianqie=sum((e_jianqie(:).^2))/(m*n);
psnr_jianqie=10*log10(255^2/mse_jianqie);%原始水印与提取水印的峰值信噪比
% disp(['攻击后峰值信噪比psnr=',num2str(psnr_jianqie)])
figure(2)
subplot(221)
imshow(cover_image);
title('原图');
subplot(222);
imshow(I0);
title('水印图');
%显示嵌入水印后的图象
subplot(223);
uint8_watermarked_image_jianqie=uint8(watermarked_image_jianqie);
imshow(uint8_watermarked_image_jianqie)
title('剪切攻击后嵌入水印图')  
subplot(224);
imshow(double(wm_jianqie));
title('剪切攻击后提取水印图')
%% 高斯噪声
I_gaosi=imnoise(I0,'gaussian',0,0.01);%高斯加噪
[wm_gaosi0,watermarked_image_gaosi,wm_gaosi]=ldpc_dct(I_gaosi,cover_image);%ldpc_dct嵌入提取
e_gaosi=wm_gaosi0-wm_gaosi;
[m,n]=size(e_gaosi);
mse_gaosi=sum((e_gaosi(:).^2))/(m*n);
psnr_gaosi=10*log10(255^2/mse_gaosi);%原始水印与提取水印的峰值信噪比
% disp(['高斯噪声攻击后峰值信噪比psnr=',num2str(psnr_gaosi)])
figure(3)
subplot(221)
imshow(cover_image);
title('原图');
subplot(222);
imshow(I0);
title('水印图');

%显示嵌入水印后的图象
subplot(223);
uint8_watermarked_image_gaosi=uint8(watermarked_image_gaosi);
imshow(uint8_watermarked_image_gaosi)
title('高斯噪声攻击后嵌入水印图')  
subplot(224);
imshow(double(wm_gaosi));
title('高斯噪声攻击后提取水印图')
%% 旋转攻击
%%9.rotate 45 旋转
I_xuanzhuan=imrotate(I0,45,'bilinear','crop');%旋转45度
[wm_xuanzhuan0,watermarked_image_xuanzhuan,wm_xuanzhuan]=ldpc_dct(I_gaosi,cover_image);%ldpc_dct嵌入提取
e_xuanzhuan=wm_xuanzhuan0-wm_xuanzhuan;
[m,n]=size(e_xuanzhuan);
mse_xuanzhuan=sum((e_xuanzhuan(:).^2))/(m*n);
psnr_xuanzhuan=10*log10(255^2/mse_xuanzhuan);%原始水印与提取水印的峰值信噪比
% disp(['旋转攻击后峰值信噪比psnr=',num2str(psnr_xuanzhuan)])
figure(4)
subplot(221)
imshow(cover_image);
title('原图');
subplot(222);
imshow(I0);
title('水印图');
function H = makeLdpc(M, N, method, noCycle, onePerCol)
% Create R = 1/2 low density parity check matrix
%
%  M        : Number of row
%  N        : Number of column
%  method   : Method for distributing non-zero element
%             {0} Evencol : For each column, place 1s uniformly at random
%             {1} Evenboth: For each column and row, place 1s uniformly at random
%  noCyle   : Length-4 cycle
%             {0} Ignore (do nothing)
%             {1} Eliminate
%  onePerCol: Number of ones per column
%
%  H        : Low density parity check matrix                   
%
%
% Copyright Bagawan S. Nugroho, 2007 
% http://bsnugroho.googlepages.com


% Number of ones per row (N/M ratio must be 2)
if N/M ~= 2
   fprintf('Code rate must be 1/2\n');
end
onePerRow = (N/M)*onePerCol;

% fprintf('Creating LDPC matrix...\n');

switch method
   % Evencol
   case {0}
      % Distribute 1s uniformly at random within column
      for i = 1:N
         onesInCol(:, i) = randperm(M)';
      end
        
      % Create non zero elements (1s) index
      r = reshape(onesInCol(1:onePerCol, :), N*onePerCol, 1);
      tmp = repmat([1:N], onePerCol, 1);
      c = reshape(tmp, N*onePerCol, 1);
      
      % Create sparse matrix H
      H = full(sparse(r, c, 1, M, N));
      
   % Evenboth
   case {1}
      % Distribute 1s uniformly at random within column
      for i = 1:N
         onesInCol(:, i) = randperm(M)';
      end
        
      % Create non zero elements (1s) index
      r = reshape(onesInCol(1:onePerCol, :), N*onePerCol, 1);
      tmp = repmat([1:N], onePerCol, 1);
      c = reshape(tmp, N*onePerCol, 1);
     
      % Make the number of 1s between rows as uniform as possible     
      
      % Order row index
      [r, ix] = sort(r);
      
      % Order column index based on row index
      for i = 1:N*onePerCol
         cSort(i, :) = c(ix(i));
      end
      
      % Create new row index with uniform weight
      tmp = repmat([1:M], onePerRow, 1);
      r = reshape(tmp, N*onePerCol, 1);
      
      % Create sparse matrix H
      % Remove any duplicate non zero elements index using logical AND
      S = and(sparse(r, cSort, 1, M, N), ones(M, N));
      H = full(S);     
      
end % switch

% Check rows that have no 1 or only have one 1
for i = 1:M
   
   n = randperm(N);
   % Add two 1s if row has no 1
   if length(find(r == i)) == 0
      H(i, n(1)) = 1;
      H(i, n(2)) = 1;
   % Add one 1 if row has only one 1   
   elseif length(find(r == i)) == 1
      H(i, n(1)) = 1;
   end

end % for i

% If desired, eliminate any length-4 cycle
if noCycle == 1
   
   for i = 1:M
      % Look for pair of row - column
      for j = (i + 1):M         
         w = and(H(i, :), H(j, :));
         c1 = find(w);
         lc = length(c1);
         if lc > 1
                       
            % If found, flip one 1 to 0 in the row with less number of 1s
            if length(find(H(i, :))) < length(find(H(j, :)))
               % Repeat the process until only one column left 
               for cc = 1:lc - 1
                  H(j, c1(cc)) = 0;
               end

  

 

  
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
  
66
  
67
  
68
  
69
  
70
  
71
  
72
  
73
  
74
  
75
  
76
  
77
  
78
  
79
  
80
  
81
  
82
  
83
  
84
  
85
  
86
  
87
  
88
  
89
  
90
  
91
  
92
  
93
  
94
  
95
  
96
  
97
  
98
  
99
  
100
  
101
  
102
  
103
  
104
  
105
  
106
  
107
  
108
  
109
  
110
  
111
  
112
  
113
  
114
  
115
  
116
  
117
  
118
  
119
  
120
  
121
  
122
  
123
  
124
  
125
  
126
  
127
  
128
  
129
  
130
  
131
  
132
  
133
  
134
  
135
  
136
  
137
  
138
  
139
  
140
  
141
  
142
  
143
  
144
  
145
  
146
  
147
  
148
  
149
  
150
  
151
  
152
  
153
  
154
  
155
  
156
  
157
  
158
  
159
  
160
  
161
  
162
  
163
  
164
  
165
  
166
  
167
  
168
  
169
  
170
  
171
  
172
  
173
  
174
  
175
  
176
  
177
  
178
  
179
  
180
  
181
  
182
  
183
  
184
  
185
  
186
  
187
  
188
  
189
  
190
  
191
  
192
  
193
  
194
  
195
  
196
  
197
  
198
  
199
  
200
  
201
  
202
  
203
  
204
  
205
  
206
  
207
 

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社，2020.
[2]杨丹,赵海滨,龙哲.MATLAB图像处理实例详解[M].清华大学出版社，2013.
[3]周品.MATLAB图像处理与图形用户界面设计[M].清华大学出版社，2013.
[4]刘成龙.精通MATLAB图像处理[M].清华大学出版社，2015.
[5]万谊丹.基于Arnold和DCT的抗剪切攻击图像水印研究[J].网络安全技术与应用. 2021,(08)

文章来源: qq912100926.blog.csdn.net，作者：海神之光，版权归原作者所有，如需转载，请联系作者。

原文链接：qq912100926.blog.csdn.net/article/details/116151831