%% Copyright (C) 2012 Telematic and Telecommunication Laboratory (LTT),       
% Dipartimento di Ingegneria dell'Informazione- Università di Siena
% via Roma 56 53100 - Siena, Italy                   
% 
%----------------------------------------------------------------------
%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
%
%----------------------------------------------------------------------
function[remapped percent] = image_remapper_recursive(img, mapping, plotFig)
%function[remapped] = image_remapper(img,mapping, plotFig)
%Takes the input image img and remaps a set of its pixels according to the
%mapping vector. Choice of remapping pairs is done basing on a perceptual
%criterion (remap high variance regions first), and update the criterion
%from time to time to lower false-boundaries artifact
%--- INPUT ARGS ---
%-> img: grayscale image
%-> mapping: 256 x 256 matrix. The i-th element of the matrix says how many
%       pixels of the i-th value are to be moved to the j-th value.
%-> plotFig [OPTIONAL] set to 'y' to have figures plotted.
%
%--- OUTPUT ARGS ---
%-> remapped: remapped image
%-> percent: percent of pixels that have been modified

%convert to grayscale if needed, warn the user
if size(img,3)>1
   img = rgb2gray(img);
   fprintf('Warning: image was not grayscale, it has been converted.\n');
end

%% Clean mapping matrix: remove redundancies
cleanCount = 0;
totCleanPix = 0;
for i=1:size(mapping,1)
    for j=1:i-1 %execute just in the lower triangle
        if mapping(i,j)>0 && mapping(j,i)>0
            if mapping(i,j)>mapping(j,i)
                mapping(i,j) = mapping(i,j) - mapping(j,i);
                totCleanPix = totCleanPix + min(mapping(j,i),mapping(i,j));
                mapping(j,i) = 0;
            else
                mapping(j,i) = mapping(j,i) - mapping(i,j);
                totCleanPix = totCleanPix + min(mapping(j,i),mapping(i,j));                
                mapping(i,j) = 0;                
            end
            cleanCount = cleanCount+1;
        end
    end
end
%fprintf('Cleaned %d remappings (%.3f%% pixels in total)',cleanCount,double(totCleanPix)/double(sum(sum(mapping))));  
%Compute map of the perceptual distortion
[varianceImg]=funcVarianceMap(double(img), 5, -Inf);

%% Apply each remapping
%
fix_img_lin = img(:);
img_lin = img(:);
var_lin = varianceImg(:)+eps;


[a srt_idxs_var] = sort(var_lin,'descend');

%Partition the sorted list in P parts
P = 4;
randomized = split_and_randomize(srt_idxs_var, P);

faults=false;
numreppix=0;
for i=1:size(mapping,1)
    srcval = i-1; %AMPL counts from 0, Matlab from 1...
    okpix = fix_img_lin(randomized)==srcval; %list available pixels for this remapping
    for j=1:size(mapping,2)    
        tbd = uint32(round(mapping(i,j)));
        if tbd>0 && j~=i
            trgval = j-1; %AMPL counts from 0, Matlab from 1...
%            fprintf('Mapping %d pixels from %d to %d\n', tbd, srcval, trgval);                                
                srcidxs = find(randomized.*okpix,tbd,'first'); %butta via gli indici che non vanno bene, trova i primi tbd 
                img_lin(randomized(srcidxs)) = trgval;  
                okpix(srcidxs)=false; %avoid re-using the same source pixels more than once.                
                replaced = length(srcidxs);
                numreppix = numreppix + replaced;
                if replaced < (tbd-1) %ignore difference of 1 pixel due to rounding of the mapping matrix
                   fprintf('Warning: image doesn''t have enough pixels of level (%d): %d out of %d remapped\n',uint8(srcval),replaced,tbd);
                   faults = true;
                   break;
                end                
        end
    end
    %update distortion criteria
    tmp = uint8(reshape(img_lin,size(img,1),size(img,2)));    
    [dummy updated] = ssim_index(tmp, img);
    updated = PadOrCrop(updated, size(img));
    [dummy randomized] = sort(updated(:),'descend');
end
remapped = uint8(reshape(img_lin,size(img,1),size(img,2)));

%Apply soft smoothing in disturbed regions
a=0;

percent = numreppix/numel(img);
if faults==false
%    fprintf('Remapping finished: %d pixels remapped (%.3f%%).\n',numreppix,100*numreppix/numel(img));
else
    fprintf('Remapping finished with some faults: %d pixels remapped (%.3f%%).\n',numreppix,100*numreppix/numel(img));
end


%% Plot if requested
if nargin>2 && strcmpi(plotFig,'y')
    psnr = psnr_mse_maxerr(img,remapped);
     try
       ssim = ssim_index(img,remapped);
     catch ME
         ssim = NaN;
         fprintf('Cannot evaluate SSIM (out of memory?)... reported error: %s.\n',ME.message);
     end
    figure();
    subplot(1,3,1);imshow(img);title('Original image');
    subplot(1,3,2);imshow(remapped);title('Remapped image');
    subplot(1,3,3);imshow(remapped-img,[]);title(sprintf('Difference image (expanded)\n. PSRN:%.3f  --  SSIM:%.3f',psnr,ssim));
    figure();
    subplot(1,2,1);imhist(img);title('Original image histogram');
    subplot(1,2,2);imhist(remapped);title('Remapped image histogram');
end

end

function[randomized] = split_and_randomize(array,P)
    psz = floor(length(array)/P);
    for j=1 : psz : length(array)-psz
        randomized(j:j+psz) = permuteRand(array(j:j+psz));
    end
    randomized = [randomized array(end-(length(array)-length(randomized))+1:end)']';
end

function[permuted] = permuteRand(array)
    assert(min(size(array))==1,'Error: cannot permutate multidemensional arrays.');
    n = length(array);
    p_idx = randperm(n);
    permuted = array(p_idx);
    %permuted = array;
end