function[d] = chi2dist(P,Q,normalize)
%% 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[d] = chi2dist(P,Q,normalize)
%Calculate Chi-square distance.
%
%--- INPUT ARGS ---
%-> P: first histogram (1-by-N array of numbers)
%-> Q: second histogram (1-by-N array of numbers, same length of P)
%-> normalize [OPTIONAL] : set to 'normalize' to normalize histograms before
%    evaluating the distance
%
%--- OUTPUT ARGS ---
%-> d: Chi-square distance between P and Q

if nargin>2 &&  strcmpi(normalize,'normalize')
    P=P/sum(P);
    Q=Q/sum(Q);
end

    if length(P)~=length(Q)
      fprintf('Error: cannot compare histograms of different length.');
      d=[];
      return
    end
    
    denom = P+Q;
    zeroidxs = denom~=0; %ignore bins that are empty in both histograms!
    d = sum( ((P(zeroidxs)-Q(zeroidxs)).^2)./(denom(zeroidxs)./2) );
    
end