function I = mutualinfo(P,m,n,i,j)

Pi = getmargX(P,m,n,i)+realmin;
Pj = getmargX(P,m,n,j)+realmin;
Pij = getmargX(P,m,n,[i j]);

I = 0;
for x=1:m
  for y=1:m
    xy = coord2ind([x y],[m m]);
    I = I + Pij(xy)*log(Pij(xy)/(Pi(x)*Pj(y)));
  end
end
return