% SUBPIX2D Sub-pixel locations in 2D image
%
% Usage: [rs, cs] = subpix2d(r, c, L);
%
% Arguments:
% r, c - row, col vectors of extrema to pixel precision.
% L - 2D corner image
%
% Returns:
% rs, cs - row, col vectors of valid extrema to sub-pixel
% precision.
%
% Note that the number of sub-pixel extrema returned can be less than the number
% of integer precision extrema supplied. Any computed sub-pixel location that
% is more than 0.5 pixels from the initial integer location is rejected. The
% reasoning is that this implies that the extrema should be centred on a
% neighbouring pixel, but this is inconsistent with the assumption that the
% input data represents extrema locations to pixel precision.
%
% The sub-pixel locations are solved by forming a Taylor series representation
% of the corner image values in the vicinity of each integer location extrema
%
% L(x) = L + dL/dx' x + 1/2 x' d2L/dx2 x
%
% x represents a position relative to the integer location of the extrema. This
% gives us a quadratic and we solve for the location where the gradient is zero
% - these are the extrema locations to sub-pixel precision
%
% Reference: Brown and Lowe "Invariant Features from Interest Point Groups"
% BMVC 2002 pp 253-262
%
% See also: SUBPIX3D
% Copyright (c) 2010 Peter Kovesi
% Centre for Exploration Targeting
% School of Earth and Environment
% The University of Western Australia
% peter.kovesi at uwa edu au
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% July 2010
function [rs, cs] = subpix2d(R, C, L)
if ndims(L) == 3
error('Corner image must be grey scale');
end
[rows, cols] = size(L);
x = zeros(2, length(R)); % Buffer for storing sub-pixel locations
m = 0; % Counter for valid extrema
for n = 1:length(R)
r = R(n); % Convenience variables
c = C(n);
% If coords are too close to boundary skip
if r < 2 || c < 2 || ...
r > rows-1 || c > cols-1
continue
end
% Compute partial derivatives via finite differences
% 1st derivatives
dLdr = (L(r+1,c) - L(r-1,c))/2;
dLdc = (L(r,c+1) - L(r,c-1))/2;
D = [dLdr; dLdc]; % Column vector of derivatives
% 2nd Derivatives
d2Ldr2 = L(r+1,c) - 2*L(r,c) + L(r-1,c);
d2Ldc2 = L(r,c+1) - 2*L(r,c) + L(r,c-1);
d2Ldrdc = (L(r+1,c+1) - L(r+1,c-1) - L(r-1,c+1) + L(r-1,c-1))/4;
% Form Hessian from 2nd derivatives
H = [d2Ldr2 d2Ldrdc
d2Ldrdc d2Ldc2 ];
% Solve for location where gradients are zero - these are the extrema
% locations to sub-pixel precision
if rcond(H) < eps
continue; % Skip to next point
% warning('Hessian is singular');
else
dx = -H\D; % dx is location relative to centre pixel
% Check solution is within 0.5 pixels of centre. A solution
% outside of this implies that the extrema should be centred on a
% neighbouring pixel, but this is inconsistent with the
% assumption that the input data represents extrema locations to
% pixel precision. Hence these points are rejected
if all(abs(dx) <= 0.5)
m = m + 1;
x(:,m) = [r;c] + dx;
end
end
end
% Extract the subpixel row and column values from x noting we just
% have m valid extrema.
rs = x(1, 1:m);
cs = x(2, 1:m);