www.gusucode.com > 个用于图像标定的算法matlab源码程序 > code16/lk20_p2/affine_ic.m

    function fit = affine_ic(img, tmplt, p_init, n_iters, var, verbose)
% AFFINE_IC - Affine image alignment using inverse-compositional algorithm
%   FIT = AFFINE_IC(IMG, TMPLT, P_INIT, N_ITERS, PERC_OCC, VERBOSE)
%   Align the template image TMPLT to an example image IMG using an
%   affine warp initialised using P_INIT. Iterate for N_ITERS
%   iterations. Parameter VAR is not used in this algorithm.
%   To display the fit graphically set VERBOSE non-zero.
%
%   p_init = [p1, p3, p5
%             p2, p4, p6];
%
%   This assumes greyscale images and rectangular templates.
%
%   c.f. Baker-Matthews

% Iain Matthews, Carnegie Mellon University, Pittsburgh


if nargin<6 verbose = 0; end
if nargin<5 error('Not enough input arguments'); end

% Common initialisation
init_a;

% Pre-computable things ---------------------------------------------------

% 3) Evaluate gradient of T
[nabla_Tx nabla_Ty] = gradient(tmplt);

% 4) Evaluate Jacobian - constant for affine warps
dW_dp = jacobian_a(w, h);

% 5) Compute steepest descent images, VT_dW_dp
VT_dW_dp = sd_images(dW_dp, nabla_Tx, nabla_Ty, N_p, h, w);

% 6) Compute Hessian and inverse
H = hessian(VT_dW_dp, N_p, w);
H_inv = inv(H);

% Algorithm -------------------------
for f=1:n_iters
  % 1) Compute warped image with current parameters
  IWxp = warp_a(img, warp_p, tmplt_pts);

  % 2) Compute error image - NB reversed
  error_img = IWxp - tmplt;
  
  % -- Save current fit parameters --
  fit(f).warp_p = warp_p;
  fit(f).rms_error = sqrt(mean(error_img(:) .^2));
  
  % -- Show fitting? --
  if verbose
    disp(['Inverse-Compositional [',num2str(f-1),']: RMS = ',num2str(fit(f).rms_error)]);
    verb_plot_a(verb_info, warp_p, tmplt_pts, error_img);
  end
  
  % -- Really iteration 1 is the zeroth, ignore final computation --
  if (f == n_iters) break; end

  % 7) Compute steepest descent parameter updates
  sd_delta_p = sd_update(VT_dW_dp, error_img, N_p, w);

  % 8) Compute gradient descent parameter updates
  delta_p = H_inv * sd_delta_p;
  
  % 9) Update warp parmaters
  warp_p = update_step(warp_p, delta_p);

  
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function warp_p = update_step(warp_p, delta_p)
% Compute and apply the update

delta_p = reshape(delta_p, 2, 3);

% Convert affine notation into usual Matrix form - NB transposed
delta_M = [delta_p; 0 0 1];	
delta_M(1,1) = delta_M(1,1) + 1;
delta_M(2,2) = delta_M(2,2) + 1;

% Invert compositional warp
delta_M = inv(delta_M);

% Current warp
warp_M = [warp_p; 0 0 1];	
warp_M(1,1) = warp_M(1,1) + 1;
warp_M(2,2) = warp_M(2,2) + 1;

% Compose
comp_M = warp_M * delta_M;	
warp_p = comp_M(1:2,:);
warp_p(1,1) = warp_p(1,1) - 1;
warp_p(2,2) = warp_p(2,2) - 1;