www.gusucode.com > 基于互信息熵的图像配准方法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 ) ) );