系列文章目录

隐写算法笔记之WOW

前言

一、失真函数的设计

 ρ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和通过小波系数改变点的像素后之间的差异进行加权：

 （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 filterhpdf = [-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 filterlpdf = (-1).^(0:numel(hpdf)-1).*fliplr(hpdf); % fliplr使矩阵X沿垂直轴左右翻转：对于1维的 hpdf，fliplr(hpdf)是hpdf的倒序% construction of 2D wavelet filtersF{1} = lpdf"*hpdf;F{2} = hpdf"*lpdf;F{3} = hpdf"*hpdf;%% Get embedding costs% inicializationcover = double(cover);p = params.p;wetCost = 10^10;sizeCover = size(cover);% add paddingpadSize = max([size(F{1})"; size(F{2})"; size(F{3})"]);coverPadded = padarray(cover, [padSize padSize], "symmetric");% compute directional residual and suitability /xi for each filterxi = 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);endrho_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 /rhomax_cost = max(cost);rho = reshape(max_cost, M_rho,N_rho);% adjust embedding costsrho(rho > wetCost) = wetCost; % threshold on the costsrho(isnan(rho)) = wetCost; % if all xi{} are zero threshold the costrhoP1 = rho;rhoM1 = rho;rhoP1(cover==255) = wetCost; % do not embed +1 if the pixel has max valuerhoM1(cover==0) = wetCost; % do not embed -1 if the pixel has min value%% Embedding simulatorstego = 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 1-eps)) = 0;        Ht = sum(H);    endendend