系列文章目录
隐写算法笔记之WOW
前言
DESIGNING STEGANOGRAPHIC DISTORTION USING DIRECTIONAL FILTERS 这篇论文是
Vojt
ech Holub和Jessica Fridrich于2012年发表的一篇会议论文。这篇论文提出了一种“累加隐写失真”的方法,即改变一个像素后,定向高通滤波器的输出变化将被加权,然后使用Holder范数进行汇总,以定义各个像素的成本。这篇笔记将简单总结论文中的一些要点。
一、失真函数的设计
ρij意为将像素Xij变为像素Yij所产生的损失。
1.方向滤波器
(1)嵌入算法应该嵌入到不容易建模的噪声区域,即不平滑的地方,例如图像的纹理区域。
(2)使用滤波器组评估多个方向的平滑度;βn由n个多方向高通滤波器组成而且核是统一的。
(3)第k个残差为,星号由卷积镜像填充。
(4)如果残差值对于个别i,j和所有k来说都很大,那么该像素在任何方向都不平滑,因此难以建模进行隐写分析。
(5)滤波器的选择:
KB:无方向滤波器。
Sobel算子:边缘检测器
WDFB-H:使用Haar的基于小波的方向滤波器
WDFB-D:使用DB-8的基于小波的方向滤波器
小波组由三个滤波器组成,K(1)K(2)K(3),可获得水平垂直对角三个方向的残差。
2.嵌入位置的选择
(1)对R和通过小波系数改变点的像素后之间的差异进行加权:
ξ称为嵌入“适用性”,是残差的绝对值与残差变化的绝对值之间的相关性;所有像素的ξ可以一次性被计算完。
(2)通过汇总所有适用性ξ来计算嵌入成本ρij。
(3)聚合规则必须满足以下两点:
a.ξ越大,ρ越小。
b.ξ为0是ρ为正无穷。
因此选取Holder范数的倒数:
(4)嵌入的更改被限制为±1。
二、实验
这里就不再描述作者的算法与其他算法的对比了,有兴趣的可以去论文中看一看,直接附上matlab的代码:
clc;clear all;
indir = 'D:\BOSSbase';
Output_path = 'D:\WOW\stego\spatial';
feature_path1 = 'D:\WOW\feature\spatial\cover_2000.mat';
% my_SRMQ1_pgm(indir,feature_path1);
params.p = -1;
payload = 0.2;
x = 0;
err = zeros(2,6);
% indir=[input '\'];
feature_path2 = ['D:\WOW\feature\spatial\stego_2000' '_WOW_' num2str(payload*100) '.mat'];
if exist(Output_path,'dir'); rmdir(Output_path,'s'); end
if ~exist(Output_path,'dir'); mkdir(Output_path); end
flist = dir([indir '\*.pgm']);
flen = length(flist);
fprintf('%s%d\n', 'the num of the files: ',flen);
tic
for weight = 0:0.1:0.5
x = x +1;
params.w = weight;
parfor i = 1:flen
fprintf('%d%s\n',i, [' processing image: ' flist(i).name]);
in_file_name = [indir '\' flist(i).name];
out_file_dir = [ Output_path '\' flist(i).name];
cover = imread(in_file_name);
[stego, dist] = WOW(cover, payload, params);
stego=uint8(stego);
imwrite(stego,out_file_dir,'pgm');
%show_s_dif(cover,stego);
end
toc
%%
my_SRMQ1_pgm(Output_path,feature_path2);
test_error = my_ensemble(feature_path1,feature_path2);
err(1,x) = payload;
err(2,x) = test_error;
end
err
save ('err_0','err');
失真函数(上一段的WOW部分):
%% Get 2D wavelet filters - Daubechies 8
% 1D high pass decomposition filter
hpdf = [-0.0544158422, 0.3128715909, -0.6756307363, 0.5853546837, 0.0158291053, -0.2840155430, -0.0004724846, 0.1287474266, 0.0173693010, -0.0440882539, ...
-0.0139810279, 0.0087460940, 0.0048703530, -0.0003917404, -0.0006754494, -0.0001174768];
% 1D low pass decomposition filter
lpdf = (-1).^(0:numel(hpdf)-1).*fliplr(hpdf); % fliplr使矩阵X沿垂直轴左右翻转:对于1维的 hpdf,fliplr(hpdf)是hpdf的倒序
% construction of 2D wavelet filters
F{1} = lpdf'*hpdf;
F{2} = hpdf'*lpdf;
F{3} = hpdf'*hpdf;
%% Get embedding costs
% inicialization
cover = double(cover);
p = params.p;
wetCost = 10^10;
sizeCover = size(cover);
% add padding
padSize = max([size(F{1})'; size(F{2})'; size(F{3})']);
coverPadded = padarray(cover, [padSize padSize], 'symmetric');
% compute directional residual and suitability \xi for each filter
xi = cell(3, 1);
Rho = cell(3, 1);
for fIndex = 1:3
% compute residual
R = conv2(coverPadded, F{fIndex}, 'same');
% show_cost_dis(R)
% compute suitability
xi{fIndex} = conv2(abs(R), rot90(abs(F{fIndex}), 2), 'same');
% correct the suitability shift if filter size is even 滤波器的大小为偶数
% 这个设计是由滤波器的size问题带来的,因为偶数的滤波器模板,确定不了中心点
if mod(size(F{fIndex}, 1), 2) == 0, xi{fIndex} = circshift(xi{fIndex}, [1, 0]); end;
if mod(size(F{fIndex}, 2), 2) == 0, xi{fIndex} = circshift(xi{fIndex}, [0, 1]); end;
% remove padding
xi{fIndex} = xi{fIndex}(((size(xi{fIndex}, 1)-sizeCover(1))/2)+1:end-((size(xi{fIndex}, 1)-sizeCover(1))/2), ((size(xi{fIndex}, 2)-sizeCover(2))/2)+1:end-((size(xi{fIndex}, 2)-sizeCover(2))/2));
Rho{fIndex} = spatial_neighbourhood(xi{fIndex}.^p);
end
rho_1 = Rho{1};
rho_2 = Rho{2};
rho_3 = Rho{3};
[M_rho,N_rho] = size(rho_1);
cost = zeros(3,M_rho*N_rho);
cost(1,:) = reshape(rho_1, 1, M_rho*N_rho);
cost(2,:) = reshape(rho_2, 1, M_rho*N_rho);
cost(3,:) = reshape(rho_3, 1, M_rho*N_rho);
% compute embedding costs \rho
max_cost = max(cost);
rho = reshape(max_cost, M_rho,N_rho);
% adjust embedding costs
rho(rho > wetCost) = wetCost; % threshold on the costs
rho(isnan(rho)) = wetCost; % if all xi{} are zero threshold the cost
rhoP1 = rho;
rhoM1 = rho;
rhoP1(cover==255) = wetCost; % do not embed +1 if the pixel has max value
rhoM1(cover==0) = wetCost; % do not embed -1 if the pixel has min value
%% Embedding simulator
stego = EmbeddingSimulator(cover, rhoP1, rhoM1, payload*numel(cover), false);
distortion_local = rho(cover~=stego);
distortion = sum(distortion_local);
%% --------------------------------------------------------------------------------------------------------------------------
% Embedding simulator simulates the embedding made by the best possible ternary coding method (it embeds on the entropy bound).
% This can be achieved in practice using "Multi-layered syndrome-trellis codes" (ML STC) that are asymptotically aproaching the bound.
function [y] = EmbeddingSimulator(x, rhoP1, rhoM1, m, fixEmbeddingChanges)
n = numel(x);
lambda = calc_lambda(rhoP1, rhoM1, m, n);
pChangeP1 = (exp(-lambda .* rhoP1))./(1 + exp(-lambda .* rhoP1) + exp(-lambda .* rhoM1));
pChangeM1 = (exp(-lambda .* rhoM1))./(1 + exp(-lambda .* rhoP1) + exp(-lambda .* rhoM1));
if fixEmbeddingChanges == 1
RandStream.setGlobalStream(RandStream('mt19937ar','seed',139187));
else
RandStream.setGlobalStream(RandStream('mt19937ar','Seed',sum(100*clock)));
end
randChange = rand(size(x));
y = x;
y(randChange < pChangeP1) = y(randChange < pChangeP1) + 1;
y(randChange >= pChangeP1 & randChange < pChangeP1+pChangeM1) = y(randChange >= pChangeP1 & randChange < pChangeP1+pChangeM1) - 1;
function lambda = calc_lambda(rhoP1, rhoM1, message_length, n)
l3 = 1e+3;
m3 = double(message_length + 1);
iterations = 0;
while m3 > message_length
l3 = l3 * 2;
pP1 = (exp(-l3 .* rhoP1))./(1 + exp(-l3 .* rhoP1) + exp(-l3 .* rhoM1));
pM1 = (exp(-l3 .* rhoM1))./(1 + exp(-l3 .* rhoP1) + exp(-l3 .* rhoM1));
m3 = ternary_entropyf(pP1, pM1);
iterations = iterations + 1;
if (iterations > 10)
lambda = l3;
return;
end
end
l1 = 0;
m1 = double(n);
lambda = 0;
alpha = double(message_length)/n;
% limit search to 30 iterations
% and require that relative payload embedded is roughly within 1/1000 of the required relative payload
while (double(m1-m3)/n > alpha/1000.0 ) && (iterations<30)
lambda = l1+(l3-l1)/2;
pP1 = (exp(-lambda .* rhoP1))./(1 + exp(-lambda .* rhoP1) + exp(-lambda .* rhoM1));
pM1 = (exp(-lambda .* rhoM1))./(1 + exp(-lambda .* rhoP1) + exp(-lambda .* rhoM1));
m2 = ternary_entropyf(pP1, pM1);
if m2 < message_length
l3 = lambda;
m3 = m2;
else
l1 = lambda;
m1 = m2;
end
iterations = iterations + 1;
end
end
function Ht = ternary_entropyf(pP1, pM1)
p0 = 1-pP1-pM1;
P = [p0(:); pP1(:); pM1(:)];
H = -((P).*log2(P));
H((P<eps) | (P > 1-eps)) = 0;
Ht = sum(H);
end
end
end
数据集下载链接:
https://dde.binghamton.edu/download/