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 &lt;= 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 &gt;= 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) &lt; 4 || length (y) &lt; 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 &lt; 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 &gt; 3 points 
i = 1;
sum = 0;
while (i &lt;= 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 &gt; 3 points 
i = 1;
sum = 0;
while (i &lt;= 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?
Just realized that this is an equation involving homogeneous vectors; thus the 3-vectors should not equal, they have the same direction but may differ in magnitude by a nonzero scale factor.