I'm doing a direct implementation of DLT algorithm but I can't figure out why, at the end the matrix H I get doesn't (even remotely) produce x2 = H*x1
The algorithm is not that complicated so I can't see what I'm missing at this stage.
%%point matching section (manual for now)
% point selection
[originalpoints1, normpoints1, translation1, scaling1 ]= pointprocessing(rootimage, false);
[originalpoints2, normpoints2, translation2, scaling2 ]= pointprocessing(sideimage, true);
if(length (originalpoints1) ~= length (originalpoints2))
error('points must be referenced in both images');
end
%% DLT transform application
% getting A
i = 1;
A = zeros(2*length(originalpoints1), 9);
while(i <= length(normpoints1))
x1 = normpoints1(1,i);
y1 = normpoints1(2,i);
w1 = normpoints1(3,i);
x2 = normpoints2(1,i);
y2 = normpoints2(2,i);
w2 = normpoints2(3,i);
a1 = - w2 * (normpoints1(:,i)');
a2 = y2 * (normpoints1(:,i)');
a3 = w2 * (normpoints1(:,i)');
a4 = - x2 * (normpoints1(:,i)');
z = zeros(1,3);
a = [z a1 a2; a3 z a4];
A(2*i-1, :) = a(1, :);
A(2*i, :) = a(2, :);
i = i + 1;
end
% getting H
[U D V] = svd(A);
h = V(:, end);
normH = reshape(h, 3, 3)';
% denormalization
T1 = scaling1 * translation1;
T2 = scaling2 * translation2;
H = inv(T2) * normH * T1;
The normalizing process works as intended - I've double checked, performs the translation and scaling (none on Z) as intended. But for clarity sake here that one also:
function [ originalpoints, normpoints, translation, scaling ] = pointprocessing( image, flag )
% point selection
figure(1);
imshow(image);
axis on;
if (~flag)
title('select n >= 4 reference points, double click last');
else
title('mark the same reference points (same order), double click last');
end
[x, y] = getpts();
close(1);
if(length (x) < 4 || length (y) < 4)
error('selected too few points');
end
% coordinate conversion - centroid calculation
center_x = mean(x);
center_y = mean(y);
originalpoints = [x'; y'; ones(1, length (x))];
% coordinate conversion - apply the translation
translation = [1 0 -center_x; 0 1 -center_y; 0 0 1];
translatedpoints = translation * originalpoints;
% coordinate conversion - apply isotropic scaling
goal = sqrt(2);
heuristic = meandistance(translatedpoints);
if(heuristic < goal)
error('average distance between points too (1 pixel) small? double check')
end
scaling_factor = (length(translatedpoints) * goal) / adaptedmeandistance(translatedpoints);
scaling = [scaling_factor 0 0; 0 scaling_factor 0; 0 0 1];
normpoints = scaling * translatedpoints;
end
function [sum] = adaptedmeandistance(point_matrix)
%assumes > 3 points
i = 1;
sum = 0;
while (i <= length(point_matrix))
sum = sum + abs(realsqrt( point_matrix(1,i)^2 + point_matrix(2,i)^2));
i = i + 1;
end
end
function [distance] = meandistance(point_matrix)
%assumes > 3 points
i = 1;
sum = 0;
while (i <= length(point_matrix))
sum = sum + euclidian(point_matrix(1,i), point_matrix(2,i), 0, 0);
i = i + 1;
end
distance = sum / length(point_matrix);
end
function [distance] = euclidian(point1_x, point1_y, point2_x, point2_y)
distance = abs(realsqrt( (point1_x - point2_x)^2 + (point1_y - point2_y)^2));
end
The issue is that I get, for instance, an H of
0.4787 0.0051 186.8102
0.0702 0.4762 144.1276
0.0000 0.0000 0.4514
Transforming originalpoints1
1.0e+03 *
0.1465 2.3465 2.3465 0.1985
0.1545 0.1345 3.0545 3.0945
0.0010 0.0010 0.0010 0.0010
into something,
1.0e+03 *
0.2577 1.3108 1.3257 0.2976
0.2280 0.3729 1.7634 1.6316
0.0005 0.0006 0.0007 0.0006
that's not originalpoints2
1.0e+03 *
0.5545 2.2865 1.9265 0.5105
0.4905 0.6505 2.5625 2.7985
0.0010 0.0010 0.0010 0.0010
The problem might be on the scale, should I manually enforce the scale of each point to 1? Doesn't that mean there's something wrong with H?
