Computing a matrix which transforms a quadrangle to another quadrangle in 2D

Question

In the figure below the goal is to compute the homography matrix H which transforms the points a1 a2 a3 a4 to their counterparts b1 b2 b3 b4. That is:

[b1 b2 b3 b4] = H * [a1 a2 a3 a4]

What way would you suggest to be the best way to calculate H(3x3). a1...b4 are points in 2D which are represented in homogeneous coordinate systems ( that is [a1_x a1_y 1]', ...). EDIT: For these types of problems we use SVD, So i would like to see how this can be simply done in Matlab.

EDIT:

Here is how I initially tried to solve it using svd (H=Q/P) in Maltlab. Cosider the following code for the given example

px=[0 1 1 0];  % a square
py=[1 1 0 0];

qx=[18 18 80 80];    % a random quadrangle
qy=[-20 20 60 -60];
if (DEBUG)
  fill(px,py,'r');
  fill(qx,qy,'r');
end

Q=[qx;qy;ones(size(qx))];
P=[px;py;ones(size(px))];
H=Q/P;
H*P-Q
answer:
   -0.0000         0         0         0         0
  -20.0000   20.0000  -20.0000   20.0000    0.0000
   -0.0000         0         0         0   -0.0000

I am expecting the answer to be a null matrix but it is not!... and that's why I asked this question in StackOverflow. Now, we all know it is a projective transformation not obviously Euclidean. However, it is good to know if in general care calculating such matrix using only 4 points is possible.

matrix computation

This is not my area, but I believe the best you can do is a least-squares solution, as you have more constraints in your equation than you have free parameters. — Isaac, Jul 30 '12 at 16:45
Thats exactly right, SVD is the solution however I am missing something somewhere in my code. — C graphics, Jul 30 '12 at 17:00
What are you trying right now? By analogy to linear least squares, my first guess would be `H = B * A' * inv( A * A' )` — Isaac, Jul 30 '12 at 17:09
in Matlab simply H=B/A where B=[b1 b2 b3 b4], A=[a1 a2 a3 a4] — C graphics, Jul 30 '12 at 17:11
And you think there should exist a matrix that does better? I'm not sure what you're asking... — Isaac, Jul 30 '12 at 17:15
Cross reference: see [this post on Math SE](http://math.stackexchange.com/a/339033/35416) for a mathematical discussion of how to find the projective transformation, without special tools from matlab. — MvG, Aug 28 '15 at 12:59

score 2 · Answer 1 · answered Jul 30 '12 at 18:24

2

You can try the cp2tform function, which infers spatial transformation from control point pairs. Since parallel is not preserved in your case, you should set the transformtype to be 'projective'. More info is here

answered Jul 30 '12 at 18:24

chaohuang

3,965
4
27
35

score 1 · Answer 2 · answered Jul 30 '12 at 17:40

1

You can use the DLT algorithm for this purpose. There are MATLAB routines available for doing that in Peter Kovesi's homepage.

answered Jul 30 '12 at 17:40

jmbr

3,298
23
23

Amro · Accepted Answer · 2012-07-31T23:43:11.883

Using the data you posted:

P = [px(:) py(:)];
Q = [qx(:) qy(:)];

Compute the transformation:

H = Q/P;

apply the transformation:

Q2 = H*P;

Compare the results:

err = Q2-Q

output:

err =
   7.1054e-15   7.1054e-15
  -3.5527e-15  -3.5527e-15
  -1.4211e-14  -2.1316e-14
   1.4211e-14   1.4211e-14

which is zeros for all intents and purposes..

EDIT:

So as you pointed out in the comments, the above method will not compute the 3x3 homography matrix. It simply solves the system of equations with as many equations as points provided:

H * A = B   -->   H = B*inv(A)   -->   H = B/A (mrdivide)

Otherwise, MATLAB has the CP2TFORM function in the image processing toolbox. Here is an example applied to the image shown:

%# read illustration image
img = imread('https://i.stack.imgur.com/ZvaZK.png');
img = imcomplement(im2bw(img));

%# split into two equal-sized images
imgs{1} = img(:,fix(1:end/2));
imgs{2} = img(:,fix(end/2:end-1));

%# extract the four corner points A and B from both images
C = cell(1,2);
for i=1:2
    %# some processing
    I = imfill(imgs{i}, 'holes');
    I = bwareaopen(imclearborder(I),200);
    I = imfilter(im2double(I), fspecial('gaussian'));

    %# find 4 corners
    C{i} = corner(I, 4);

    %# sort corners in a consistent way (counter-clockwise)
    idx = convhull(C{i}(:,1), C{i}(:,2));
    C{i} = C{i}(idx(1:end-1),:);
end

%# show the two images with the detected corners
figure
for i=1:2
    subplot(1,2,i), imshow(imgs{i})
    line(C{i}(:,1), C{i}(:,2), 'Color','r', 'Marker','*', 'LineStyle','none')
    text(C{i}(:,1), C{i}(:,2), num2str((1:4)'), 'Color','r', ...
        'FontSize',18, 'Horiz','left', 'Vert','bottom')
end

With the corners detected, now we can obtain the spatial transformation:

%# two sets of points
[A,B] = deal(C{:});

%# infer projective transformation using CP2TFORM
T = cp2tform(A, B, 'projective');

%# 3x3 Homography matrix
H = T.tdata.T;
Hinv = T.tdata.Tinv;

%# align points in A into B
X = tformfwd(T, A(:,1), A(:,2));

%# show result of transformation
line(X([1:end 1],1), X([1:end 1],2), 'Color','g', 'LineWidth',2)

The result:

>> H = T.tdata.T
H =
      0.74311    -0.055998    0.0062438
      0.44989      -1.0567   -0.0035331
      -293.31       62.704      -1.1742

>> Hinv = T.tdata.Tinv
Hinv =
       -1.924     -0.42859   -0.0089411
      -2.0585      -1.2615   -0.0071501
       370.68       39.695            1

screenshot

We can confirm the calculation ourselves:

%# points must be in Homogenous coordinates (x,y,w)
>> Z = [A ones(size(A,1),1)] * H;
>> Z = bsxfun(@rdivide, Z, Z(:,end))   %# divide by w
Z =
          152           57            1
          219          191            1
           62          240            1
           92          109            1

which maps to the points in B:

%# maximum error
>> max(max( abs(Z(:,1:2)-B) ))
ans =
   8.5265e-14

Thats not what I was looking for. In your solution H is 4x4, which means if we incorporate n points, H will become nxn. We are looking for a 3x3 Homography matrix.- but Thanks anyways — C graphics, Jul 31 '12 at 20:24
@Cgraphics: that's true, thanks for pointing that out. Perhaps since only 4 points are enough to solve the system, it shouldn't be a problem... Either way, I have an example using [CP2TFORM](http://www.mathworks.com/help/toolbox/images/ref/cp2tform.html) function as [chaohuang](http://stackoverflow.com/a/11727305/97160) suggested. It will correctly infer the 3x3 Homography matrix of the projective transformation from a minimum of 4 control points. I'll post it shortly, maybe someone will find it useful — Amro, Jul 31 '12 at 23:02

score 0 · Answer 4 · answered Jul 30 '12 at 19:00

0

Combine all the coordinates into column vectors. In 2D case they are: ax, ay, bx, by;

Define:

A = [ax, ay];
B = [bx, by];
input = [ones(size(A, 1), 1), A];

Calculate transformation matrix:

H = input \ B;

If you need to apply it on new observations A_new:

input_new = [ones(size(A_new, 1), 1), A_new];
B_new = input_new * H;

Edit: you can also calculate H by: H = inv(input' * input) * input' * B; but \ does if safely.

answered Jul 30 '12 at 19:00

Serg

13,470
8
36
47

@Amro: Sorry to hijack this thread, but I wanted Amro's opinion. If you look at my answer, I have referenced my solution to a previous question. I believe it's not _technically_ an exact duplicate but perhaps the difference is irrelevant in SO. What do you think? – Jacob Jul 30 '12 at 20:33
@Jacob: absolutely, that's a great answer. Perhaps you should post code here adapted to this problem. (btw, there seem to be a dead link pointing to some slides). IMO these questions are similar although not exact duplicates. – Amro Jul 30 '12 at 21:17

score 0 · Answer 5 · edited May 23 '17 at 12:19

0

I discussed a related problem in an answer to another SO question (includes a MATLAB solution). In the linked problem, the output quadrilateral was a rectangle, but my answer handles the general problem.

edited May 23 '17 at 12:19

Community

1
1

answered Jul 30 '12 at 20:29

Jacob

34,255
14
110
165

1

Thanks Jacob, You are right it is exactly same question. I was looking for Homography in fact. – C graphics Jul 30 '12 at 22:59

Computing a matrix which transforms a quadrangle to another quadrangle in 2D

5 Answers5

EDIT:

Linked