www.pudn.com > 1111111registrationbasedoncorrelation.rar > GMI.m

function G = GMI(Is, Ir, sigma) 
% % % -- referrence by  
% % %   "Image registration by maximization of combined mutual information 
% % %       and gradient information" -- J P W Pluim, ... 
% % %  
% % % -- cal gradient of (Is, Ir) 
% % %  
% % % -- algrithm 
% % % 1. DetXs, DetXr = conv image Is and Ir with first derivatives of Gaussian kernel 
% % % 2. AlphaXsr = the angle between the DetXs and DetXr 
% % % 3. WAlpha = the weight for angle 
% % % 4. G = sum( WAlpha * min(|DetXs|, |DetXr|) ) 
 
% % ============== main begins ==================== 
% % % 1. DetXs, DetXr = conv image Is and Ir with first derivatives of Gaussian kernel 
if nargin <3 
    sigma =0.5; 
end 
% sigma = 1/3; 
% siz = [3 3]; 
% [DetXsx DetXsy MagXs] = calDetX(Is, sigma, siz); 
% [DetXrx DetXry MagXr] = calDetX(Ir, sigma, siz); 
[DetXsx DetXsy MagXs] = convgau(Is, sigma); 
[DetXrx DetXry MagXr] = convgau(Ir, sigma); 
 
 
% % % 2. AlphaXsr = the angle between the DetXs and DetXr 
AlphaXsr = acos( ( DetXsx .* DetXrx + DetXsy .* DetXry ) ... 
    ./( MagXs .* MagXr +0.00001 ) ); 
 
% % % 3. AlphaW = the weight for angle 
AlphaW = ( cos(2* AlphaXsr) +1 )/2; 
 
% % % 4. G = sum( AlphaW * min(|DetXs|, |DetXr|) ) 
 
G = AlphaW * min(MagXs, MagXr)' ; 
% G= sum(G(:)); 
 
% % ============== main ends ==================== 
 
% % -- sub functions begins 
% % -- 
function [ax,ay,mag]= convgau(I, sigma) 
a = im2double(I); 
  % Magic numbers 
  GaussianDieOff = .0001;   
  PercentOfPixelsNotEdges = .7; % Used for selecting thresholds 
  ThresholdRatio = .4;          % Low thresh is this fraction of the high. 
   
  % Design the filters - a gaussian and its derivative 
   
  pw = 1:30; % possible widths 
  ssq = sigma*sigma; 
  width = max(find(exp(-(pw.*pw)/(2*sigma*sigma))>GaussianDieOff)); 
  if isempty(width) 
    width = 1;  % the user entered a really small sigma 
  end 
 
  t = (-width:width); 
  gau = exp(-(t.*t)/(2*ssq))/(2*pi*ssq);     % the gaussian 1D filter 
 
% arg   = -(t.*t)/(2*ssq); 
% gau   = exp(arg); 
 
  % Find the directional derivative of 2D Gaussian (along X-axis) 
  % Since the result is symmetric along X, we can get the derivative along 
  % Y-axis simply by transposing the result for X direction. 
  [x,y]=meshgrid(-width:width,-width:width); 
  %dgau2D=-x.*exp(-(x.*x+y.*y)/(2*ssq))/(pi*ssq); 
dgau2D = - x.* exp(-(x.*x+y.*y)/(2*ssq)) /ssq; 
 
  % Convolve the filters with the image in each direction 
  % The canny edge detector first requires convolution with 
  % 2D gaussian, and then with the derivitave of a gaussian. 
  % Since gaussian filter is separable, for smoothing, we can use  
  % two 1D convolutions in order to achieve the effect of convolving 
  % with 2D Gaussian.  We convolve along rows and then columns. 
 
%   %smooth the image out 
%   aSmooth=imfilter(a,gau,'conv','replicate');   % run the filter accross rows 
%   aSmooth=imfilter(aSmooth,gau','conv','replicate'); % and then accross columns 
%    
  %apply directional derivatives 
aSmooth = a;   
  ax = imfilter(aSmooth, dgau2D, 'conv','replicate'); 
  ay = imfilter(aSmooth, dgau2D', 'conv','replicate');%转制还是需要的！！ 
% X = [ax ay]; 
ax= ax(:)'; ay =ay(:)'; 
  mag = sqrt((ax.*ax) + (ay.*ay)); 
   
%   magmax = max(mag(:)); 
%   if magmax>0 
%     mag = mag / magmax;   % normalize 
%   end 
%    
  mag = mag(:)';