stc/exp/steganalysers/aump/aump.m

   1 function beta = aump(X,m,d)
   2 %
   3 % AUMP LSB detector as described by L. Fillatre, "Adaptive Steganalysis of
   4 % Least Significant Bit Replacement in Grayscale Natural Images", IEEE
   5 % TSP, October 2011.
   6 %
   7 % X = image to be analyzed
   8 % m = pixel block size
   9 % d = q - 1 = polynomial degree for fitting (predictor)
  10 % beta = \hat{\Lambda}^\star(X) detection statistic
  11 %
  12
  13 X = double(X);
  14 [Xpred,~,w] = Pred_aump(X,m,d);       % Polynomial prediction, w = weights
  15 r = X - Xpred;                        % Residual
  16 Xbar = X + 1 - 2 * mod(X,2);          % Flip all LSBs
  17 beta = sum(sum(w.*(X-Xbar).*r));      % Detection statistic
  18
  19
  20 function [Xpred,S,w] = Pred_aump(X,m,d)
  21 %
  22 % Pixel predictor by fitting local polynomial of degree d = q - 1 to
  23 % m pixels, m must divide the number of pixels in the row.
  24 % OUTPUT: predicted image Xpred, local variances S, weights w.
  25 %
  26 % Implemention follows the description in: L. Fillantre, "Adaptive
  27 % Steganalysis of Least Significant Bit Replacement in Grayscale Images",
  28 % IEEE Trans. on Signal Processing, 2011.
  29 %
  30
  31 sig_th = 1;               % Threshold for sigma for numerical stability
  32 q = d + 1;                % q = number of parameters per block
  33 Kn = numel(X)/m;          % Number of blocks of m pixels
  34 Y = zeros(m,Kn);          % Y will hold block pixel values as columns
  35 S = zeros(size(X));       % Pixel variance
  36 Xpred = zeros(size(X));   % Predicted image
  37
  38 H = zeros(m,q);           % H = Vandermonde matrix for the LSQ fit
  39 x1 = (1:m)/m;
  40 for i = 1 : q, H(:,i) = (x1').^(i-1); end
  41
  42 for i = 1 : m             % Form Kn blocks of m pixels (row-wise) as
  43     aux = X(:,i:m:end);   % columns of Y
  44     Y(i,:) = aux(:);
  45 end
  46
  47 p = H\Y;                  % Polynomial fit
  48 Ypred = H*p;              % Predicted Y
  49
  50 for i = 1 : m             % Predicted X
  51     Xpred(:,i:m:end) = reshape(Ypred(i,:),size(X(:,i:m:end))); % Xpred = l_k in the paper
  52 end
  53
  54 sig2 = sum((Y - Ypred).^2) / (m-q);           % sigma_k_hat in the paper (variance in kth block)
  55 sig2 = max(sig_th^2 * ones(size(sig2)),sig2); % Assuring numerical stability
  56 % le01 = find(sig2 < sig_th^2);
  57 % sig2(le01) = (0.1 + sqrt(sig2(le01))).^2;   % An alternative way of "scaling" to guarantee num. stability
  58
  59 Sy = ones(m,1) * sig2;                        % Variance of all pixels (order as in Y)
  60
  61 for i = 1 : m             % Reshaping the variance Sy to size of X
  62     S(:,i:m:end) = reshape(Sy(i,:),size(X(:,i:m:end)));
  63 end
  64
  65 s_n2 = Kn / sum(1./sig2);                     % Global variance sigma_n_bar_hat^2 in the paper
  66 w = sqrt( s_n2 / (Kn * (m-q)) ) ./ S;         % Weights
  67