Jose Delpiano

Andrew, many thanks. I'm sorry about my delay. I added the improvements you suggested. The code is much faster now.

Divyesh Varade, mutual information is a measure based on pixel statistics, so it can't be applied pixel by pixel. Applying it to small windows could be a good idea for some applications.

Great Job ! However, I'm just wondering why you're using hist in mi function and histc in hist2 function. I would have used histc everywhere. Is there any difference, or it is an arbitrary choice ?

Thanks for the code it was very helpful. However, I want to compare two images using MI. If I do it pixel by pixel using the given code, I mostly get a warning: divide by zero and receive a value NaN for MI. Now I can try and work with very small windows also but that would seem to reduce the quality of my results. Is there another way or by enhancing this code somehow to analyse the two images together almost pixelwise.

the current implementation of both MI and hist2 are quite slow;

function n=fast_hist2(A,B,L)
ma=min(A(:));
MA=max(A(:));
mb=min(B(:));
MB=max(B(:));

% For sensorimotor variables, in [-pi,pi]
% ma=-pi;
% MA=pi;
% mb=-pi;
% MB=pi;

% Scale and round to fit in {0,...,L-1}
A=round((A-ma)*(L-1)/(MA-ma+eps));
B=round((B-mb)*(L-1)/(MB-mb+eps));
n=zeros(L);
x=0:L-1;
for i=0:L-1
    n(i+1,:) = histc(B(A==i),x,1);
end
end

is about 10x quicker, and gives the same result (assuming that the input is sensible, i believe you get different behavior if you pass in all nan's, or anything which would cause matlabs 'hist' function to error) ; given that the mi function calls hist on both the inputs individually first, this should not be a problem.

also there is need for a bugfix in mi, 256 bins are hard coded, irrespective of what you pass in.

the minf function also does unnecessary allocation.

Using the functions above and below rather than the attached the running time is about 22 seconds rather than 240, for my set of MI calculations. (2500 calls for MI between two 50x50 windows)

function I=fast_mi(A,B,L)
%MI Determines the mutual information of two images or signals
%
% I=mi(A,B)
%
% Assumption: 0*log(0)=0
%
% See also WENTROPY.

% jfd, 15-11-2006
% 01-09-2009, added case of non-double images

A=double(A);
B=double(B);
    
na = fast_hist(A(:),L);
na = na/sum(na);

nb = fast_hist(B(:),L);
nb = nb/sum(nb);

n2 = fast_hist2(A,B,L);
n2 = n2/sum(n2(:));

I=sum(minf(n2,na'*nb));
%I=sum(I);

% -----------------------

function y=minf(pab,papb)

I=find(papb(:)>1e-12 & pab(:)>1e-12); % function support
y=pab(I).*log2(pab(I)./papb(I));