function[remapped] = histogram_universal_antiforensics(img, maxDisc)
%% Copyright (C) 2012 Telematic and Telecommunication Laboratory (LTT),       
% Dipartimento di Ingegneria dell'Informazione- Università di Siena
% via Roma 56 53100 - Siena, Italy                   
% 
% Code authors: M. Barni, M. Fontani, B. Tondi
%----------------------------------------------------------------------
%Permission to use, copy, or modify this software and its documentation
%for educational and research purposes only and without fee is hereby
%granted, provided that this copyright notice and the original authors'
%names appear on all copies and supporting documentation. This program
%shall not be used, rewritten, or adapted as the basis of a commercial
%software or hardware product without first obtaining permission of the
%authors. The authors make no representations about the suitability of
%this software for any purpose. It is provided "as is" without express
%or implied warranty.
%----------------------------------------------------------------------
%
%This is an implementation of the algorithm for removing statistical
%traces left into the histogram of an image by any processing operator. 
%Please refer to the following paper:
%
%M. Barni, M. Fontani, B. Tondi, "A Universal Technique to Hide Traces of 
%Histogram-Based Image Manipulations", ACM Int. Workshop on Multimedia and 
%Security (MMSEC) 2012, Coventry, UK, 2012.
%
%Kindly report any suggestions or corrections to marco.fontani@unisi.it
%
%----------------------------------------------------------------------
%Notice: this code requires AMPL and ipopt to solve the optimization
%problem. AMPL must be in the system path.
%
%--- INPUT ARGS ---
%-> img: image (N-by-M matrix of 8-bit integers) that should be 
%       processed (notice: this code works only for grayscale images. 
%       Image will be converted to grayscale from RGB if needed).
%-> maxDisc: maximum displacement allowed (typical values 
%       range from 2 to 8), refer to the paper.
%
%--- OUTPUT ARGS ---
%-> remapped: image resulting from the application of
%       counter-forensic algorithm.
%
%Usage example:
%   remapped = histogram_universal_antiforensics(myImage, 4);
%
%See the results:
%   imshow(remapped)  %Shows the produced image.
%
%========================================================================


    DBfile = 'MIRFLICKR.mat'; %file containing histograms of untouched images (use generateDB.m to create)
    addpath(genpath('./'));
    
        %% **************** STEP 0 *****************
        % *             Prepare image              *
        % ******************************************
        %convert to grayscale if needed
        if size(img,3)>1
            img = rgb2gray(img);
            warning('\nImage is RGB, converting to grayscale automatically...\n');%#ok<WNTAG>
        end
        hist_img = hist(img(:),0:255);        

        %% **************** STEP 1 *****************
        % *            Search in DB               *
        % *****************************************
        %Find the 10 best matching histograms among those in the DB
        fprintf('Searching best match in the histogram DB... ')
        [matching_hist] = MMfind_best_match(img, 10, DBfile, '');
        fprintf('\t[DONE]\n');
        
        fprintf('Starting optimization with Dmax = %d... ',maxDisc);
        %% **************** STEP 2 *****************
        % *         Find optimal mapping          *
        % *****************************************
        %Choose the mapping that yields best (minimum) obj. fnc. value.
        fncval=[];mapps=cell(1,1);Nimg=[];Normaliz=[];
        for i=1:size(matching_hist,1)
            %Prepare optimization: generate optimization problem data (AMPL)
            Nimg(i) = sum(hist_img);
            Normaliz(i) = generateFileDat_1(hist_img, matching_hist{i}, maxDisc);
            %Use AMPL to write problem and run the solver
            [s retstr] = system('ampl Divergence.run');
            %fprintf('SYSTEM MESSAGE WHEN RUNNING AMPL:\n\t%s\n',retstr);
            %Read results
            [mapps{i} fncval(i)]= readFileDat_1('n_opt.dat');     %#ok<*AGROW>
        end
        [dummy IDXbest] = min(fncval);
        mapping = mapps{IDXbest};
        %rescale mapping from empiric probabilities to actual pixel values
        mapping = uint32((mapping./Normaliz(IDXbest)).*Nimg(IDXbest));
        hist_target = matching_hist{IDXbest};
        fprintf('\t[DONE]\n');
        
        
        %% **************** STEP 3 *****************
        % *            Apply mapping              *
        % *****************************************
        fprintf('Starting pixel remapping (this may take a while)...');        
        [remapped] = image_remapper_recursive(img, mapping, 'n');
        %get histogram of remapped image
        hist_remapped = hist(remapped(:),0:1:255);
        fprintf('\t[DONE]\n');
        
        %% **************** STEP 4 *****************
        % *      Check quality of output image     *
        % *****************************************
        %evaluate PSNR ans SSIM between remapped and attacked images
        psnr = psnr_mse_maxerr(img,remapped);
        try
            ssim = ssim_index(img,remapped);
        catch ME
            ssim = NaN;
            fprintf('[Warning]: Cannot evaluate SSIM (out of memory?)... reported error: %s.\n',ME.message);
        end
        %% **************** PLOT and CHATS *****************
        % *      Plot histograms and similar..     *
        % *****************************************
        %Plot histograms
        figure();set(gcf,'OuterPosition',[0 400 1200 400]);
        subplot(1,3,1); bar(hist_img);title('Input histogram');
        subplot(1,3,2); bar(hist_target);title('Target (from DB)');
        subplot(1,3,3); bar(hist_remapped);title('Output histogram');
        
        %Plot images
        figure();set(gcf,'OuterPosition',[0 0 1200 400])
        subplot(1,3,1);imshow(uint8(img));title('Input image');
        subplot(1,3,2);imshow(remapped);title('Output image');
        subplot(1,3,3);imshow(remapped-uint8(img),[]);title(sprintf('Difference image (expanded)\n. PSRN:%.3f  --  SSIM:%.3f',psnr,ssim));
        
        %Print log to screen
        fprintf('**** Image successfully remapped! ****\nSimilarity:\t PSNR = %.3f\tSSIM = %.3f\n', psnr, ssim);

end
