function [nmi,jointEntropy] = mutualInfo(ima,imb)

% MUTUALINFO calculates the normalized Mutual Information (MI) of two 
% images using their histograms.
% 
% 
% Inputs:
%           ima     :  First image
%           imb     :  Second image
% 
% Outputs:
%           nmi             :  Normalized mutual information
%           jointEntropy	:  Joint entropy
% 
% 
% (C)	Mohammad Haghighat, University of Miami
%       haghighat@ieee.org


if nargin < 2   % Check correct number of arguments
    error('There should be two input images!');
end

if size(ima) ~= size(imb)   % Check if the input images are of the same size
    error('Size of the input images must be the same!');
end

[m,n] = size(ima);
ima = double(ima);
imb = double(imb);

% Change a, b to the range of [1 max] (maximum gray-level of the image)
if min(ima(:)) < 1
    ima = ima - min(ima(:)) + 1;
end
if min(imb(:))<1
    imb = imb - min(imb(:)) + 1;
end

% Calculate the maximum gray-levels of the images
aMax = max(ima(:));
bMax = max(imb(:));

% Calculate image histograms
aHist = zeros(1,aMax);
bHist = zeros(1,bMax);
jointHist = zeros(aMax,bMax);

for i=1:m
    for j=1:n
        
        indexa =  ima(i,j);
        indexb =  imb(i,j);
                
        aHist(indexa) = aHist(indexa) + 1;
        bHist(indexb) = bHist(indexb) + 1;
        
        jointHist(indexa,indexb) = jointHist(indexa,indexb) + 1;
        
   end
end

% Calculate entropies
histSum = sum(sum(aHist));
index = find(aHist~=0);
p = aHist/histSum;  % Probability density function (pdf) of a
aEntropy = sum(sum(-p(index).*log2(p(index))));	% Entropy of a

histSum = sum(sum(bHist));
index = find(bHist~=0);
p = bHist/histSum;  % pdf of b
bEntropy = sum(sum(-p(index).*log2(p(index))));	% Entropy of b

histSum = sum(sum(jointHist));
index = find(jointHist~=0);
p = jointHist/histSum;  % Joint pdf of a & b
jointEntropy = sum(sum(-p(index).*log2(p(index))));	% Joint Entropy between a & b

% Calculate the mutual information (mi) of a & b
mi = aEntropy + bEntropy - jointEntropy;	% Mutual Info.
nmi = 2 * mi ./ (aEntropy + bEntropy);      % Normalized mutual info.

end
