41
$\begingroup$

I went through this thread Mapping Irregular Quadrilateral to a Rectangle

If i know the 4 corresponding points in image say
p1->p1'
p2->p2'
p3->p3'
p4->p4'

then how to compute pi(x,y) from pi'(x,y)

enter image description here enter image description here

i don't know how to compute elements in Homography matrix H from those 8 known points
[x']= [h11 h12 h13] [x]
[y']= [h21 h22 h23] [y]
[(1)]=[h31 h32 (1)] [(1)]

[Excuse me. I am not sure if I should extend this question, or create a new one, since I can't post comments on threads]

I want to ask the same question, but using absolute values so I can visualize it. Lets say my points on the image plane are:

p[0] = x:407 y:253
p[1] = x:386 y:253
p[2] = x:406 y:232
p[3] = x:385 y:232

these points are in a 500px width x 333px height image plane with 0,0 at top left corner. These points represents a picture of a real plane where a 30mm side square are located. Assuming this picture was taken by a fixed camera at origin heading Z axis.

So, I know the physical distance between p0,p1 ; p0,p2 ; p1,p3; p2,p3 are 30mm.

But is it possible to get the X,Y,Z from each of these points using only this information above?

$\endgroup$
1
  • $\begingroup$ Thank you guyz. after searching a lot i found this wiimote librarry by brian brianpeek.com/page/wiimotelib. in his code warp matrix calculation is implemented. that's what i was actually looking for. $\endgroup$
    – Sagar
    Mar 23, 2014 at 19:28

4 Answers 4

35
$\begingroup$

You can compute the homography matrix H with your eight points with a matrix system such that the four correspondance points $(p_1, p_1'), (p_2, p_2'), (p_3, p_3'), (p_4, p_4')$ are written as $2\times9$ matrices such as:

$p_i = \begin{bmatrix} -x_i \quad -y_i \quad -1 \quad 0 \quad 0 \quad 0 \quad x_ix_i' \quad y_ix_i' \quad x_i' \\ 0 \quad 0 \quad 0 \quad -x_i \quad -y_i \quad -1 \quad x_iy_i' \quad y_iy_i' \quad y_i' \end{bmatrix}$

It is then possible to stack them into a matrix $P$ to compute:

$PH = 0$

Such as:

$PH = \begin{bmatrix} -x_1 \quad -y_1 \quad -1 \quad 0 \quad 0 \quad 0 \quad x_1x_1' \quad y_1x_1' \quad x_1' \\ 0 \quad 0 \quad 0 \quad -x_1 \quad -y_1 \quad -1 \quad x_1y_1' \quad y_1y_1' \quad y_1' \\ -x_2 \quad -y_2 \quad -1 \quad 0 \quad 0 \quad 0 \quad x_2x_2' \quad y_2x_2' \quad x_2' \\ 0 \quad 0 \quad 0 \quad -x_2 \quad -y_2 \quad -1 \quad x_2y_2' \quad y_2y_2' \quad y_2' \\ -x_3 \quad -y_3 \quad -1 \quad 0 \quad 0 \quad 0 \quad x_3x_3' \quad y_3x_3' \quad x_3' \\ 0 \quad 0 \quad 0 \quad -x_3 \quad -y_3 \quad -1 \quad x_3y_3' \quad y_3y_3' \quad y_3' \\ -x_4 \quad -y_4 \quad -1 \quad 0 \quad 0 \quad 0 \quad x_4x_4' \quad y_4x_4' \quad x_4' \\ 0 \quad 0 \quad 0 \quad -x_4 \quad -y_4 \quad -1 \quad x_4y_4' \quad y_4y_4' \quad y_4' \\ \end{bmatrix} \begin{bmatrix}h1 \\ h2 \\ h3 \\ h4 \\ h5 \\ h6 \\ h7 \\ h8 \\h9 \end{bmatrix} = 0$

While adding an extra constraint $|H|=1$ to avoid the obvious solution of $H$ being all zeros. It is easy to use SVD $P = USV^\top$ and select the last singular vector of $V$ as the solution to $H$.

Note that this gives you a DLT (direct linear transform) homography that minimizes algebraic error. This error is not geometrically meaningful and so the computed homography may not be as good as you expect. One typically applies nonlinear least squares with a better cost function (e.g. symmetric transfer error) to improve the homography.

Once you have your homography matrix $H$, you can compute the projected coordinates of any point $p(x, y)$ such as:

$\begin{bmatrix} x' / \lambda \\ y' / \lambda \\ \lambda \end{bmatrix} = \begin{bmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{bmatrix}. \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}$

$\endgroup$
8
  • $\begingroup$ The 8x2 p_i matrices seems to be 9x2 in the answer. $\endgroup$
    – anorm
    Oct 14, 2015 at 7:29
  • $\begingroup$ @anorm They should be 9x2. I think 8x2 is just a typo $\endgroup$
    – marcman
    Nov 9, 2015 at 19:20
  • 2
    $\begingroup$ @marcman h33 has been set to 1. That why its 8. $\endgroup$ Nov 16, 2015 at 10:31
  • 6
    $\begingroup$ @tony can you clarify the dimension mismatch, and how h1-h8 relate to h11-h33 ? $\endgroup$
    – Guig
    Jun 19, 2017 at 21:19
  • 1
    $\begingroup$ @Cypher This is an application of linear algebra. For purely learning linear algebra, I found Linear Algebra and Its Applications by David C, Lay a good resource. As a side note, this comment is off-topic to this question and should be taken offline. $\endgroup$
    – Miket25
    Aug 13, 2020 at 22:37
8
$\begingroup$

Rather than taking SVD of the system, there is an easier way to solve this. since the last component in the $H$ vector is considered as $1$, we can add one more constraint equation $(h9=1)$ (since $h_{33}$ is considered $1$) to this set and convert $P$ matrix into $9\times 9$ matrix and system will become non-homogenous.

$PH = \begin{bmatrix} -x_1 & -y_1 & -1 & 0 & 0 & 0 & x_1x_1' & y_1x_1' & x_1' \\ 0 & 0 & 0 & -x_1 & -y_1 & -1 & x_1y_1' & y_1y_1' & y_1' \\ -x_2 & -y_2 & -1 & 0 & 0 & 0 & x_2x_2' & y_2x_2' & x_2' \\ 0 & 0 & 0 & -x_2 & -y_2 & -1 & x_2y_2' & y_2y_2' & y_2' \\ -x_3 & -y_3 & -1 & 0 & 0 & 0 & x_3x_3' & y_3x_3' & x_3' \\ 0 & 0 & 0 & -x_3 & -y_3 & -1 & x_3y_3' & y_3y_3' & y_3' \\ -x_4 & -y_4 & -1 & 0 & 0 & 0 & x_4x_4' & y_4x_4' & x_4' \\ 0 & 0 & 0 & -x_4 & -y_4 & -1 & x_4y_4' & y_4y_4' & y_4' \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\ \end{bmatrix} \begin{bmatrix}h1 \\ h2 \\ h3 \\ h4 \\ h5 \\ h6 \\ h7 \\ h8 \\h9 \end{bmatrix} = \begin{bmatrix}0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\1 \end{bmatrix}$

This is a $Ax=b$ kind of problem which can be solved much faster than computing SVD.

$\endgroup$
2
  • 2
    $\begingroup$ There are some special cases where the ninth element is zero $\endgroup$ May 10, 2019 at 2:57
  • 1
    $\begingroup$ Yes, $h9$ can get close to 0. A much better way is to assume $h1 + h2 + ... + h9 = 1$. $\endgroup$ Dec 14, 2021 at 15:49
3
$\begingroup$

After a while I am answering my own question (in a way I can understand.I hope it can help other people too)

I am really sorry for not having a good math basis, but there is a GAP between information most people provide from copy/pasted formulas found on google and what I can understand. I see, by the way, many formulas have some mistakes that confuses those in the search of an answer.

So,

Basically, you have 2 collections of 4 points, A and B rectangles.

A = { p1x,p1y ; p2x,p2y ; p3x,p3y ; p4x,p4y }

B = { p1x,p1y ; p2x,p2y ; p3x,p3y ; p4x,p4y }

I want to map a point in rectangle A to rectangle B.

I then build a 8 x n matrix (M1), where n is number of points * 2(x,y). (result will be a 8x8 matrix)

for each point:
{xA, yA, 1, 0, 0, 0, -xA*xB, -yA*xB}
{0, 0, 0, xA, yA, 1, -xA*yB, -yA*yB}

I also make a 1 x 8 matrix (M2) with the "target"(B) points. Result must be 1x8 and not 8x1.

for each point:
{xB}
{yB}

Next, we are finding the so called H matrix. I am using this formula:

H = ( M1 transpose * M1 )inv * ( M1 transpose * M2)

Result is a 3 x 3 matrix (H).

Finally, any 2d point in rectangle A can be found in rectangle B using this operation:

point_in_A = (x,y,1)
tempMatrix (1x3) = H(3x3) * point_in_A(1x3)
tempMatrix will be (x, y, scale);

using tempMatrix values below:

result xy_in_B = (x / scale , y / scale);

That is it.

$\endgroup$
6
  • 11
    $\begingroup$ Please, use Latex. $\endgroup$ Jul 15, 2016 at 7:35
  • 4
    $\begingroup$ Hi, this is unfortunately not the best solution. You're finding the 3x3 homography as a solution to over-specified linear system in eight unknowns, assuming the ninth is one (also missing in your answer). This works as long as the ninth element of H is nonzero. But there are special cases when it is zero and then this method fails. Better use the SVD trick for total least squares (refer to e.g. the Hartley-Zisserman book for algorithm and details). $\endgroup$
    – the swine
    Jun 20, 2017 at 13:05
  • 3
    $\begingroup$ @MithleshUpadhyay OP may have a Latex allergy. In this case, SE's MathJax support may be a good alternate :-) $\endgroup$
    – uhoh
    Sep 12, 2017 at 0:32
  • 2
    $\begingroup$ the end result is H is a 1x8 matrix not a 3x3 $\endgroup$ Nov 3, 2017 at 20:25
  • $\begingroup$ Just to confirm that this works when actually programming it $\endgroup$
    – XLVII
    Nov 1, 2022 at 13:00
0
$\begingroup$

Explaining how to obtain the linear system of equations:

From $P^{\prime}=HP$, we have

$$\begin{bmatrix}u_i\\ v_i\\ 1\end{bmatrix}=\begin{bmatrix}h_1 & h_2 & h_3\\ h_4 & h_5 & h_6\\ h_7 & h_8 & h_9\end{bmatrix}\begin{bmatrix}x_i\\ y_i\\ 1\end{bmatrix}$$

from which we get

$\begin{align} u_i=h_1x_i+h_2y_i+h_3 & \quad (1)\\ v_i=h_4x_i+h_5y_i+h_6 & \quad (2)\\ 1=h_7x_i+h_8y_i+h_9 & \quad (3)\\ \end{align}$

  • $(1) / (3)$ and rearranging gives $h_7x_iu_i+h_8y_iu_i+h_9u_i-h_1x_i-h_2y_i-h_3=0$
  • $(2) / (3)$ and rearranging gives $h_7x_iv_i+h_8y_iv_i+h_9v_i-h_4x_i-h_5y_i-h_6=0$

Hence, each pair of points $p_i^{\prime} = \begin{bmatrix}u_i\\ v_i\\ 1\end{bmatrix}, p_i = \begin{bmatrix}x_i\\ y_i\\ 1\end{bmatrix}$, we get $2$ equations.

From $4$ pair of points we get linear system with $2\times 4 = 8$ equations, for $8$ unknowns $h_1,h_2,\ldots,h_8$.

Projective Transform (Homography) having $8$ degrees of freedom, $h_9$ is known or can be obtained from the unknown coefficients, so we have to solve for $8$ unknowns.

If we assume $h_9=1$, then we have the following linear system of equation ($A\boldsymbol{h}=\boldsymbol{b}$) to solve:

enter image description here

If we assume $\sum\limits_{i=1}^{9}h_i=1$, then we have the following linear system of equation

enter image description here

Assuming $A$ is full-rank, the system has a unique solution and can be obtained as $\boldsymbol{h}=A^{-1}\boldsymbol{b}$, or, more generally, as $\boldsymbol{h}=A^{\dagger}\boldsymbol{b}$, using the Moore-Penrose Pseudoinverse.

Now, let's demonstrate with python:

import numpy as np
import matplotlib.pylab as plt

src = np.array([[14,28],
               [288,28],
                [14,218],
               [288,218]])
dst = np.array([[34,207],
               [281,130],
                [16,29],
               [263,25]])

x, y, u, v = src[:,0], src[:,1], dst[:,0], dst[:,1]
A = np.zeros((9,9))
j = 0
for i in range(4):
    A[j,:] = np.array([-x[i], -y[i], -1, 0, 0, 0, x[i]*u[i], y[i]*u[i], u[i]])
    A[j+1,:] = np.array([0, 0, 0, -x[i], -y[i], -1, x[i]*v[i], y[i]*v[i], v[i]])
    j += 2
A[8, 8] = 1   # assuming h_9 = 1
b = [0]*8 + [1]

H = np.reshape(np.linalg.solve(A, b), (3,3))
print(H)

# [[  1.68471   ,  -0.09453488,  14.58392177],
   [  0.05006887,  -0.97074686, 242.75163122],
   [  0.00264368,   0.00027783,   1.        ]]

ps, p1s = [], [] 

for i in range(4):
    p, q = src[i], dst[i]
    p_hom = [p[0], p[1],1] 
    p1_hom = H@p_hom
    p1 = [x / p1_hom[2] for x in p1_hom[:-1]]
    ps.append(p)
    p1s.append(p1)
    print("P:{}, P':{}, H@P:{}".format(p, q, p1))

# P:[14 28], P':[ 34 207], H@P:[34.0, 207.00000000000003]
# P:[288  28], P':[281 130], H@P:[281.0, 130.0]
# P:[ 14 218], P':[16 29], H@P:[16.0, 29.0]
# P:[288 218], P':[263  25], H@P:[263.0, 24.999999999999986]

# plot the original and transformed points
ps, p1s = np.array(ps), np.array(p1s)
max_y = 300
plt.figure(figsize=(10,3))
plt.subplot(121),   plt.scatter(ps[:,0], max_y-ps[:,1], color='g', s=50), plt.title('P', size=20), plt.grid()
plt.subplot(122),   plt.scatter(p1s[:,0], max_y-p1s[:,1], color='r', s=50), plt.title(r'P$^{\prime}$', size=20), plt.grid()
plt.show()

enter image description here

As can be seen, the transformed coordinates are same as the ground-truth ones (upto numerical precision) with the estimated homography matrix.

Reference: https://www.cs.cmu.edu/~16385/lectures/lecture9.pdf

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .