<p>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. </p>

CC BY-SA 3.0
2016-08-09T18:08:38.317Z
2016-08-09T18:08:38.317Z
Kronephon
<p>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</p>

<p>The algorithm is not that complicated so I can't see what I'm missing at this stage.</p>

<pre><code>%%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;
</code></pre>

<p>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:</p>

<pre><code>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
</code></pre>

<p>The issue is that I get, for instance, an H of </p>

<pre><code>0.4787    0.0051  186.8102
0.0702    0.4762  144.1276
0.0000    0.0000    0.4514
</code></pre>

<p>Transforming originalpoints1</p>

<pre><code>   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
</code></pre>

<p>into something,</p>

<pre><code>   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
</code></pre>

<p>that's not originalpoints2</p>

<pre><code>   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
</code></pre>

<p>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?</p>

<image-processing><matlab>
DLT algorithm (4+ points) not working, what's missing?
Kronephon
# DLT algorithm (4+ points) not working, what's missing?

(Asked by: username_0 on 2016-08-09)

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?

**Tags:** image-processing, matlab


---

## 1 answers

### Answer by username_0 on 2016-08-09. Score: 3 | Accepted answer

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.


---

