www.gusucode.com > 二维重建三维重建源码程序 > 二维重建三维重建源码程序/demons3d/findupdate.m
%% Find update between two images % Changed: Dec 31st, 2011 % function [ux,uy,uz] = findupdate(F,M,vx,vy,vz,sigma_i,sigma_x) % Get Transformation [sx,sy,sz] = expfield(vx,vy,vz); % Interpolate updated image M_prime = iminterpolate(M,sx,sy,sz); % intensities at updated points area = size(M,1)*size(M,2)*size(M,3); % area of moving image % image difference diff = F - M_prime; % fixed image gradient [gx_f,gy_f,gz_f] = gradient(F); % image gradient normg2_f = gx_f.^2 + gy_f.^2 + gz_f.^2; % squared norm of gradient % moving image gradient [gx,gy,gz] = gradient(M_prime); % image gradient normg2 = gx.^2 + gy.^2 + gz.^2; % squared norm of gradient % update is Idiff / (||J||^2+(Idiff^2)/sigma_x^2) J, with Idiff = F(x)-M(x+s), and J = Grad(M(x+s)); scale = diff ./ (normg2 + diff.^2*sigma_i^2/sigma_x^2); scale(normg2==0) = 0; scale(diff ==0) = 0; scale = scale .* ((scale>=0) + (scale<0) .* sign(gx.*gx_f + gy.*gy_f + gz.*gz_f)); % avoid collapsing gradients (change sign if moving goes backward) ux = gx .* scale; uy = gy .* scale; uz = gz .* scale; % Zero non overlapping areas %ux(F==0) = 0; uy(F==0) = 0; %ux(M_prime==0) = 0; uy(M_prime==0) = 0; end