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

    function results = test_weighted(tdata, pt_offsets, alg_list, n_iters, n_tests, n_freq_tests, ...
				spatial_sigma, image_pixel_sigma, tmplt_pixel_sigma, max_spatial_error, verbose)
% TEST_WEIGHTED - Test algorithms using weighted L2 norm

% tdata has two fields:
%   tdata.img       unit8 greyscale image
%   tdata.tmplt     [x_start y_start x_end y_end] template rectangle corners
%                   matlab coordinates (1, 1) is top-left, start < end.
%
% pt_offsets is a N x 6 matrix of (x1, y1, x2, y2, x3, y3) deltas for each
% of the three affine control points. Zero mean, unit variance, 
% e.g. pt_offsets = randn(N, 6);
%
% Noise is added to the image with linearly changing std as specified in image_pixel_sigma
% Noise is added to the template with linearly changing std as specified in tmplt_pixel_sigma
%
% alg_list is a cellstr list of algorithms to run, e.g.:
% alg_list = {'affine_fa' 'affine_fc' 'affine_ia' 'affine_ic'};
%
% The "template" is created by cutting out a distorted version of tdata.tmplt.
% The "target" is the image defined by tdata.tmplt.

% Iain Matthews, Ralph Gross, Simon Baker
% Carnegie Mellon University, Pittsburgh

if nargin<11 error('Not enough input arguments'); end

% Target corner points (i.e. correct answer)
target_pts = [tdata.tmplt(1), tdata.tmplt(2);
              tdata.tmplt(1), tdata.tmplt(4);
	      tdata.tmplt(3), tdata.tmplt(4);
	      tdata.tmplt(3), tdata.tmplt(2)]';

% Target affine warp control points - triangle on the rectangle
target_affine = [tdata.tmplt(1), tdata.tmplt(2);
                 tdata.tmplt(3), tdata.tmplt(2);
		 tdata.tmplt(1) + ((tdata.tmplt(3) - tdata.tmplt(1)) / 2) - 0.5, tdata.tmplt(4)]';

% Template image dimensions
template_nx = tdata.tmplt(3) - tdata.tmplt(1) + 1;
template_ny = tdata.tmplt(4) - tdata.tmplt(2) + 1;

% Template corner points (unperturbed, rectangle)
template_pts = [1, 1;
		1, template_ny;
		template_nx, template_ny;
		template_nx, 1]';

% Template affine warp control points
template_affine = [1, 1;
		   template_nx, 1;
		   template_nx / 2, template_ny]';

% Initial warp parameters. Unperturbed translation
p_init = zeros(2, 3);
p_init(1, 3) = tdata.tmplt(1) - 1;
p_init(2, 3) = tdata.tmplt(2) - 1;

% Translate by 0.5 pixels to avoid identity warp. Warping introduces a little
% smoothing and this avoids the case where the first iteration for a forwards
% algorithm is on the "unsmoothed" unwarped image
p_init(1, 3) = p_init(1, 3) + 0.5;
p_init(2, 3) = p_init(2, 3) + 0.5;

% Scale point offsets to have required sigma
pt_offsets = pt_offsets * spatial_sigma;

% Need image to be doubles
tdata.img = double(tdata.img);

% Space for results
results = [];

% Test counters
go = 1; offset_idx = 1; all_alg_converged = 0;

% Convergence counters in field 1
for l=1:length(alg_list)
  results = setfield(results, {1}, alg_list{l}, {1}, 'n_diverge', 0);
  results = setfield(results, {1}, alg_list{l}, {1}, 'n_converge', 0);
end

% Index of successfull convergence tests
results.all_converged_idx = [];

% Might not be doing any convergence tests at all
if n_tests > 0
  cnvg_testing = 1;
else
  cnvg_testing = 0;
  cnvg_result = [];
end

% Run
while go
  if cnvg_testing
    disp(['Convergence Test: ', num2str(all_alg_converged + 1), ' (of total ', num2str(offset_idx), ')']);
  else
    disp(['Divergence Test: ', num2str(offset_idx)]);
  end 
  
  % Test points: apply current point offset to target points
  test_pts = target_affine + reshape(pt_offsets(offset_idx,:), 2, 3);
  
  % Solve for affine warp
  M = [template_affine; ones(1,size(template_affine,2))]' \ [test_pts; ones(1,size(test_pts,2))]';
  M = M';
  
  % Warp original image to get test "template" image
  tmplt_img = quadtobox(tdata.img, template_pts, M, 'bilinear');
  
  %Template w/o noise
  tmplt_clean = tmplt_img;
  
  % Add noise to template
  if ~isempty(tmplt_pixel_sigma)
    lowerB = tmplt_pixel_sigma(1);
    upperB = tmplt_pixel_sigma(2);
    sigVec = [lowerB:(upperB-lowerB)/(size(tmplt_img,2)-1):upperB];
    
    for colC=1:size(tmplt_img,2)
      tmplt_img(:,colC) = tmplt_img(:,colC)+randn(size(tmplt_img,1),1)*sigVec(colC);
    end
  end
  
  % Add noise to image
  if ~isempty(image_pixel_sigma)
    
    tmpXStart = tdata.tmplt(1);
    tmpXEnd   = tdata.tmplt(3);
    tmpYStart = tdata.tmplt(2);
    tmpYEnd   = tdata.tmplt(4);
    tmpHeight = tmpYEnd-tmpYStart+1;
    
    lowerB = image_pixel_sigma(1);
    upperB = image_pixel_sigma(2);
    sigVec = [lowerB:(upperB-lowerB)/(tmpXEnd-tmpXStart):upperB];

    noisy_img = tdata.img;

    for colC=0:(tmpXEnd-tmpXStart)
      noisy_img(tmpYStart:tmpYEnd,tmpXStart+colC) = noisy_img(tmpYStart:tmpYEnd,tmpXStart+colC)+randn(tmpHeight,1)*sigVec(colC+1);
    end
  else
    noisy_img = tdata.img;
  end
  
  var.tmplt_pixel_sigma = tmplt_pixel_sigma;
  % Initial error in affine points. This is not quite sqrt(mean(pt_offset(offset_idx,:) .^ 2)) due to p_init
  rms_pt_init = ComputePointError(test_pts, template_affine, p_init);
  
  % Run each algorithm
  for l=1:length(alg_list)
    string = ['tic; fit = ', alg_list{l}, '(noisy_img, tmplt_img, p_init, n_iters, var, verbose); t = toc;'];
    eval(string);
    
    % Evaluate point spatial error for each iteration for convergence tests
    if cnvg_testing
      rms_pt_err = zeros(length(fit), 1);
      for f=1:length(fit)
	rms_pt_error(f) = ComputePointError(test_pts, template_affine, fit(f).warp_p);
      end
      
      % Only need final spatial rms for divergence tests. Don't need fitting results
    else
      rms_pt_error = ComputePointError(test_pts, template_affine, fit(end).warp_p);
      fit = [];
    end	
    
    % Save spatial errors
    results = setfield(results, {1}, alg_list{l}, {offset_idx}, 'rms_pt_error', rms_pt_error);
    
    % Save full fitting results
    results = setfield(results, {1}, alg_list{l}, {offset_idx}, 'fit', fit);
    
    % Save fitting time
    results = setfield(results, {1}, alg_list{l}, {offset_idx}, 'time', t);
  end
  
  % Evaluate final spatial errors for all algorithms
  disp('----------------------------------------------------');
  disp(['Initial spatial rms = ',num2str(rms_pt_init)]);
  
  all_cnvg_check = 1;
  
  for l=1:length(alg_list)
    final_rms_pt_error = getfield(results, {1}, alg_list{l}, {offset_idx}, 'rms_pt_error');
    final_rms_pt_error = final_rms_pt_error(end);

    string = [alg_list{l}, ': final spatial rms = ', num2str(final_rms_pt_error)];
    
    % Convergence test
    if final_rms_pt_error > max_spatial_error
      string = [string, '   **DIVERGED**'];
      
      % Update divergence counter in first field
      n_diverge = getfield(results, {1}, alg_list{l}, {1}, 'n_diverge');
      n_diverge = n_diverge + 1;
      results = setfield(results, {1}, alg_list{l}, {1}, 'n_diverge', n_diverge);
      
      % One or more algorithms diverged
      all_cnvg_check = 0;
      
    else
      % Update convergence counter in first field
      n_converge = getfield(results, {1}, alg_list{l}, {1}, 'n_converge');
      n_converge = n_converge + 1;
      results = setfield(results, {1}, alg_list{l}, {1}, 'n_converge', n_converge);
    end
    
    disp(string);
  end
  disp('----------------------------------------------------');
  
  % Everything converged?
  if all_cnvg_check
    all_alg_converged = all_alg_converged + 1;
    
    % Update index of fully converged test results
    results.all_converged_idx = [results.all_converged_idx; offset_idx];
  end
  
  % Finished convergence testing?
  if all_alg_converged >= n_tests
    cnvg_testing = 0;
    
    % Extra tests for divergence?
    if n_freq_tests > n_tests
      if offset_idx >= n_freq_tests
	go = 0;
      end
      
      % Or just stop
    else
      go = 0;
    end
  end
  
  % Always increment index into point offsets
  offset_idx = offset_idx + 1;
  if offset_idx > size(pt_offsets, 1)
    disp('Ran out of tests!');
    return;
  end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rms_pt_error = ComputePointError(test_pts, template_affine, warp_p)
% Compute affine points rms error

% Affine for this iteration
M = [warp_p; 0 0 1];
M(1,1) = M(1,1) + 1; 
M(2,2) = M(2,2) + 1;

% Affine points
iteration_pts = M * [template_affine; ones(1, 3)];

% Error in affine points
diff_pts = test_pts - iteration_pts(1:2,:);
diff_pts = diff_pts(:);
rms_pt_error = sqrt(mean(diff_pts .^ 2));