GMI.m
上传用户:xueli1969
上传日期:2022-07-27
资源大小:19k
文件大小:3k
- 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(:)';
-