function f = extractHCFCom(image)
%EXTRACTHCFCOM
%
%   Extracts the center of mass of the adjacency histogram characteristic
%   function for steganalysis purpose.
%
%   f = extractHCFCom(image)
%
%   Input:
%     - image: a 2 dimensional 8-bit grayscale image (colour images are not
%     supported but a for loop will do the trick)
%
%   Output:
%     - f: scalar feature value
%
%   Reference paper:
%   A. D. Ker, "Steganalysis of LSB matching in grayscale images", in IEEE
%   Signal Processing Letters, vol. 12, no. 6, pp. 441-444, June 2005.
%
%   Implementation by Giacomo Cancelli, 2008
%                                                          

    % Initialize the output
    f = [];

    % Get the image dimensions
    [r c d] = size(image);
    % Checking if image type is supported
    if d ~= 1,
        disp('ERROR: You have been told that colour images were not supported!');
        return;
    end
    % Checking that image dimensions are even
    if (mod(r,2)==1 || mod(c,2)==1)
        disp('ERROR: Image dimensions should be even');
        return;
    end

    % Cast the image in double
    image = double(image);

    % Check that the value of the pixels are between 0 and 255
    % If you are using non 8-bit images you have to change the HCF2 function
    if (max(image(:))>255 || min(image(:))<0)
        disp('ERROR: The input image should be and 8 bit image.');
        return;
    end

    % Prepare the resized image for calibration
    [r c] = size(image);
    cal = image(1:2:r,1:2:c) + ...
          image(1:2:r,2:2:c) + ...
          image(2:2:r,1:2:c) + ...
          image(2:2:r,2:2:c);
    cal = floor(cal/4);

    % Compute the histogram characteristic function of the adjacency histogram
    % of the original image
    h_orig = HCF2(image);
    % Compute the centre of mass of the histogram chracteristic function
    c_orig = COM_2(h_orig);

    % Compute the histogram characteristic function of the adjacency histogram
    % of the calibrated image
    h_cal = HCF2(cal);
    % Compute the centre of mass of the histogram chracteristic function
    c_cal = COM_2(h_cal);

    % Return the ratio of the two centres of mass
    f = c_orig / c_cal;
end

%% Compute the adjacency histogram characteristic function of an image
function H = HCF2(im)
    % Compute the horizontal adjacency histogram
    h = im*256 + circshift(im,[0 -1]);
    h = histc(h(:)',[0:256*256-1]);
    h = reshape(h,[256,256])';

    % Apply a 2-dimensional Fourier transform
    HC = fft2(h);
    
    % Retrieve the norm of one quadrant
    H = abs(HC(1:size(HC,1)/2,1:ceil((size(HC,2)-1)/2)+1));
end

%% Compute the center of mass of an adjacency matrix
function c = COM_2(matrix)
    % Generate an auxiliary matrix whose elements are the sum of the row
    % columns index
    [r c] = size(matrix);
    aux = ones(r,1)*[1:c] + (ones(c,1)*[1:r])' - 2;
    % Compute the COM as defined in the paper
    C = sum(sum(aux.*matrix)) / sum(sum(matrix));
end