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(:)';