MATLAB 互信息函数

function [mi,mi1]=CalculateMI(a,b,gray_level,normalization)%,a_nogray,b_nogray) %{

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

参数：

a --输入图像a b --输入图像b

mi --图像a 和b互信息值

normalization --是否归一化互信息值，默认0不归一化 功能：

利用一个联合直方图，计算两幅图像的互信息值

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %}

% if nargin< 6

% b_nogray = []; % end

% if nargin< 5

% a_nogray = []; % end

if nargin< 4 || isempty(normalization) normalization = 0; end

if nargin< 3 || isempty(gray_level) gray_level = 256; end

a=double(a); b=double(b);

if max(max(a)) > 255 %保证图像为8位灰度图 a = a/max(max(a))*255; end

if max(max(b)) > 255

b = b/max(max(b))*255; end

if gray_level ~= 256%当灰度等级不为256时 a = a/max(max(a))*gray_level; b = b/max(max(b))*gray_level; end

[M,N]=size(a);

hab=zeros(gray_level,gray_level); if max(max(a))~=min(min(a))

a=(a-min(min(a)))/(max(max(a))-min(min(a))); else

a=zeros(M,N);

end

if max(max(b))~=min(min(b))

b=(b-min(min(b)))/(max(max(b))-min(min(b))); else

b=zeros(M,N); end

a=double(int16(a*(gray_level-1)))+1; b=double(int16(b*(gray_level-1)))+1;

for i=1:M for j=1:N index_x=a(i,j); index_y=b(i,j);

hab(index_x,index_y)=hab(index_x,index_y)+1; end end

habsum=M*N;

%%%%%%%%%%%%%%%%%%%%%%%% % if ~isempty(a_nogray)

% hab(1:a_nogray+1,:) = 0; % end

% if ~isempty(b_nogray)

% hab(:,b_nogray+1) = 0; % end

%%%%%%%%%%%%%%%%%%%%%%%% index=find(hab~=0); pab=hab/habsum;

Hab=sum(sum(-pab(index).*log2(pab(index))));

pa=sum((pab')); index=find(pa~=0);

Ha=sum(-pa(index).*log2(pa(index)));

pb=sum(pab);

index=find(pb~=0);

Hb=sum(-pb(index).*log2(pb(index)));

if normalization == 1 mi = (Ha+Hb)/Hab; else

mi=Ha+Hb-Hab;

end

mi1 = (Ha+Hb)/Hab;