# DLT algorithm (4+ points) not working, what's missing?

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

%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?