GMI.m
上传用户:xueli1969
上传日期:2022-07-27
资源大小:19k
文件大小:3k
源码类别:

通讯编程文档

开发平台:

Matlab

  1. function G = GMI(Is, Ir, sigma)
  2. % % % -- referrence by 
  3. % % %   "Image registration by maximization of combined mutual information
  4. % % %       and gradient information" -- J P W Pluim, ...
  5. % % % 
  6. % % % -- cal gradient of (Is, Ir)
  7. % % % 
  8. % % % -- algrithm
  9. % % % 1. DetXs, DetXr = conv image Is and Ir with first derivatives of Gaussian kernel
  10. % % % 2. AlphaXsr = the angle between the DetXs and DetXr
  11. % % % 3. WAlpha = the weight for angle
  12. % % % 4. G = sum( WAlpha * min(|DetXs|, |DetXr|) )
  13. % % ============== main begins ====================
  14. % % % 1. DetXs, DetXr = conv image Is and Ir with first derivatives of Gaussian kernel
  15. if nargin <3
  16.     sigma =0.5;
  17. end
  18. % sigma = 1/3;
  19. % siz = [3 3];
  20. % [DetXsx DetXsy MagXs] = calDetX(Is, sigma, siz);
  21. % [DetXrx DetXry MagXr] = calDetX(Ir, sigma, siz);
  22. [DetXsx DetXsy MagXs] = convgau(Is, sigma);
  23. [DetXrx DetXry MagXr] = convgau(Ir, sigma);
  24. % % % 2. AlphaXsr = the angle between the DetXs and DetXr
  25. AlphaXsr = acos( ( DetXsx .* DetXrx + DetXsy .* DetXry ) ...
  26.     ./( MagXs .* MagXr +0.00001 ) );
  27. % % % 3. AlphaW = the weight for angle
  28. AlphaW = ( cos(2* AlphaXsr) +1 )/2;
  29. % % % 4. G = sum( AlphaW * min(|DetXs|, |DetXr|) )
  30. G = AlphaW * min(MagXs, MagXr)' ;
  31. % G= sum(G(:));
  32. % % ============== main ends ====================
  33. % % -- sub functions begins
  34. % % --
  35. function [ax,ay,mag]= convgau(I, sigma)
  36. a = im2double(I);
  37.   % Magic numbers
  38.   GaussianDieOff = .0001;  
  39.   PercentOfPixelsNotEdges = .7; % Used for selecting thresholds
  40.   ThresholdRatio = .4;          % Low thresh is this fraction of the high.
  41.   
  42.   % Design the filters - a gaussian and its derivative
  43.   
  44.   pw = 1:30; % possible widths
  45.   ssq = sigma*sigma;
  46.   width = max(find(exp(-(pw.*pw)/(2*sigma*sigma))>GaussianDieOff));
  47.   if isempty(width)
  48.     width = 1;  % the user entered a really small sigma
  49.   end
  50.   t = (-width:width);
  51.   gau = exp(-(t.*t)/(2*ssq))/(2*pi*ssq);     % the gaussian 1D filter
  52. % arg   = -(t.*t)/(2*ssq);
  53. % gau   = exp(arg);
  54.   % Find the directional derivative of 2D Gaussian (along X-axis)
  55.   % Since the result is symmetric along X, we can get the derivative along
  56.   % Y-axis simply by transposing the result for X direction.
  57.   [x,y]=meshgrid(-width:width,-width:width);
  58.   %dgau2D=-x.*exp(-(x.*x+y.*y)/(2*ssq))/(pi*ssq);
  59. dgau2D = - x.* exp(-(x.*x+y.*y)/(2*ssq)) /ssq;
  60.   % Convolve the filters with the image in each direction
  61.   % The canny edge detector first requires convolution with
  62.   % 2D gaussian, and then with the derivitave of a gaussian.
  63.   % Since gaussian filter is separable, for smoothing, we can use 
  64.   % two 1D convolutions in order to achieve the effect of convolving
  65.   % with 2D Gaussian.  We convolve along rows and then columns.
  66. %   %smooth the image out
  67. %   aSmooth=imfilter(a,gau,'conv','replicate');   % run the filter accross rows
  68. %   aSmooth=imfilter(aSmooth,gau','conv','replicate'); % and then accross columns
  69. %   
  70.   %apply directional derivatives
  71. aSmooth = a;  
  72.   ax = imfilter(aSmooth, dgau2D, 'conv','replicate');
  73.   ay = imfilter(aSmooth, dgau2D', 'conv','replicate');%转制还是需要的!!
  74. % X = [ax ay];
  75. ax= ax(:)'; ay =ay(:)';
  76.   mag = sqrt((ax.*ax) + (ay.*ay));
  77.   
  78. %   magmax = max(mag(:));
  79. %   if magmax>0
  80. %     mag = mag / magmax;   % normalize
  81. %   end
  82. %   
  83.   mag = mag(:)';
  84.