通讯编程文档

开发平台：
Matlab

RMI.m：源码内容
							function rmi = RMI(A, B, d)
% % -- referrence:
% %       "Image similarity using mutual information of regions", D B Russakoff
% % -- Algorithm
% % Given A and B, which are m*n images
% %     1. create pi=vij for [Aij, Bij], then N = (m-2r)*(n-2r) is represented by P=[p1..pN]
% %     2. P0 = P - mean(P);
% %     3. C = P0*P0'/N, covariance matrix
% %     4. Hg(C) = log( (2pie)^(d/2) * det(C)^(d/2) )
% %     5. Hg(CA), Hg(CB), where CA is d/2*d/2 top left matrix of C, CB is bottom right
% %     6. RMI = Hg(CA)+Hg(CB)- Hg(C)
% % 
% %  -- Our method
% %     differrent in region, use only p[i+1,j] and p[i,j+1] for p[i,j]
% % 
% % % % test begins
% load t2.mat;
% % Is0 , It0
% [rt ct]=size(It0);
% [rs cs] = size(Is0);
% if rt>rs || ct>cs
%     x1 = (ct-cs)/2+1;
%     x2 = cs+x1-1; 
%     It0=It0([x1:x2],[x1:x2]);
% end
% A = Is0; 
% B = It0;
% % % % test ends
% % ======================== main begins ======================
[m n] = size(A);
% % === 1. create pi=vij for [Aij, Bij], then N = (m-2r)*(n-2r) is represented by P=[p1..pN]
if nargin <3
    d = 9;  %default set
end
P = crP(A, B, d);
N = size(P,2);
% % === 2. P0 = P - mean(P);
P0 = double(P) - repmat( mean(P,2), 1, N );    %sum(P(:))/N
% % === 3. C = P0*P0'/N, covariance matrix
C = P0*P0'/N;
% % === 4. Hg(C) = log( (2pie)^(d/2) * det(C)^(d/2) )
HgC = calHg(C, d*2);
% % === 5. Hg(CA), Hg(CB), where CA is d/2*d/2 top left matrix of C, CB is bottom right
CA = C(1:d, 1:d);
[m1 n1] = size(C);
CB = C((m1-d+1):m1, (n1-d+1):m1);
HgCA = calHg(CA, d);
HgCB = calHg(CB, d);
% % === 6. RMI = Hg(CA)+Hg(CB)- Hg(C)
rmi = HgCA + HgCB - HgC;
% % ======================== main ends ===================================
% % === subfunc1. create pi=vij for [Aij, Bij], then N = (m-2r)*(n-2r) is represented by P=[p1..pN]
function P = crP(A, B, d)
P = [];[m, n] = size(A);
% % fast algrithm
if d==3
    i = [1:(m-1)]; j = [1:(n-1)];
    Pa1 = A(i,   j)';
    Pa2 = A(i+1, j)';
    Pa3 = A(i,   j+1)';
    Pb1 = B(i,   j)';
    Pb2 = B(i+1, j)';
    Pb3 = B(i,   j+1)';
    P = [ Pa1(:)'; Pa2(:)'; Pa3(:)';...
          Pb1(:)'; Pb2(:)';Pb3(:)'];
end
if d==4
    i = [1:(m-1)]; j = [1:(n-1)];
    Pa1 = A(i,   j)';
    Pa2 = A(i+1, j)';
    Pa3 = A(i,   j+1)';
    Pa4 = A(i+1, j+1)';
    
    Pb1 = B(i,   j)';
    Pb2 = B(i+1, j)';
    Pb3 = B(i,   j+1)';
    Pb4 = B(i+1, j+1)';
%     
%     Pa1 = A([1:(m-1)], [1:(n-1)])';
%     Pa2 = A([2: m],    [1:(n-1)])';
%     Pa3 = A([1:(m-1)], [2: n]   )';
%     Pa4 = A([2: m],    [2: n]   )';
%     Pb1 = B([1:(m-1)], [1:(n-1)])';
%     Pb2 = B([2: m],    [1:(n-1)])';
%     Pb3 = B([1:(m-1)], [2: n]   )';
%     Pb4 = B([2: m],    [2: n]   )';
    P = [ Pa1(:)'; Pa2(:)';Pa3(:)'; Pa4(:)';...
          Pb1(:)'; Pb2(:)';Pb3(:)'; Pb4(:)'];
end
if d==9
    i = [2:(m-1)]; j = [2:(n-1)];
    
    Pa1 = A(i-1, j-1)';
    Pa2 = A(i,   j-1)';
    Pa3 = A(i+1, j-1)';
    Pa4 = A(i-1, j  )';
    Pa5 = A(i,   j   )';
    Pa6 = A(i+1, j  )';
    Pa7 = A(i-1, j+1)';
    Pa8 = A(i,   j+1)';
    Pa9 = A(i+1, j+1)';
        
    Pb1 = B(i-1, j-1)';
    Pb2 = B(i,   j-1)';
    Pb3 = B(i+1, j-1)';
    Pb4 = B(i-1, j  )';
    Pb5 = B(i,   j   )';
    Pb6 = B(i+1, j  )';
    Pb7 = B(i-1, j+1)';
    Pb8 = B(i,   j+1)';
    Pb9 = B(i+1, j+1)';
    P = [ Pa1(:)'; Pa2(:)';Pa3(:)'; Pa4(:)';...
          Pa5(:)'; Pa6(:)';Pa7(:)'; Pa8(:)'; Pa9(:)';...
          Pb1(:)'; Pb2(:)';Pb3(:)'; Pb4(:)';
          Pb5(:)'; Pb6(:)';Pb7(:)'; Pb8(:)'; Pb9(:)'];
end
function HgC = calHg(C, d)
det( C )
(2*pi*exp(1))^(d/2)  *sqrt( det( C ) )
HgC = log( (2*pi*exp(1))^(d/2)  *sqrt( det( C ) ) );

						