www.pudn.com > 1111111registrationbasedoncorrelation.rar > 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 ) ) );