系列文章目錄

隱寫演算法筆記之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和通過小波系數改變點的像素后之間的差異進行加權：

因此選取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