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

    function fit = affine_ic_weighted(img, tmplt, p_init, n_iters, var, verbose)
% AFFINE_IC_WEIGHTED - Affine image alignment using the inverse-compositional algorithm with a weighted norm
%   FIT = AFFINE_IC(IMG, TMPLT, P_INIT, N_ITERS, VAR, VERBOSE)
%   Align the template image TMPLT to an example image IMG using an
%   affine warp initialised using P_INIT. Iterate for N_ITERS iterations.
%   VAR.TMPLT_PIXEL_SIGMA contains the lower and upper bound of the spatially varying
%   noise assumed to be present in the template. To display the fit graphically set VERBOSE non-zero.
%
%   p_init = [p1, p3, p5
%             p2, p4, p6];
%
%   This assumes greyscale images and rectangular templates.


% Iain Matthews, Ralph Gross, Simon Baker
% 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);

% 3b) Compute weight
weight = compute_weight (size(tmplt), var.tmplt_pixel_sigma);

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

% 5a) Compute *weighted* steepest descent images, VT_dW_dp_weight
VT_dW_dp_weight = sd_images_weight(dW_dp, nabla_Tx, nabla_Ty, weight, N_p, h, w);

% 5b) 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_weight(VT_dW_dp, weight, 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_weight, 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;