www.gusucode.com > 配准结果检验程序 gobad.m对图像进行“平移”和“旋转” > peizhun/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(:)';