Calculate Homography with and without SVD

I've rendered an example for this question. With the trick of https://math.stackexchange.com/a/2619023/741822 I was able to calculate what seems to be the homography matrix (it passed 2 tests on $$p_5$$ and $$p_6$$).

Unfortunately, I'm not able to get the result via SVD which is why I created that question

3D Data

$$p_0 = (0|0|0)$$

$$p_1 = (-7|0|0)$$

$$p_2 = (3|-6|0)$$

$$p_3 = (7|-4|0)$$

$$p_4 = (3|2|0)$$

$$p_5 = (6|1|0)$$

$$p_6 = (4|7|0)$$

2D Data

$$p_0 = (700|288)$$

$$p_1 = (93|63)$$

$$p_2 = (293|868)$$

$$p_3 = (1207|998)$$

$$p_4 = (1218|309)$$

$$p_5 = (1540|502)$$

$$p_6 = (1679|128)$$

Calculation

General Formula

$$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$$

Add new row to get a $$9 \times9$$ matrix

$$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}$$

Fill in values

$$P = \displaystyle \left[\begin{matrix}-93 & -63 & -1 & 0 & 0 & 0 & -651 & -441 & -7\\0 & 0 & 0 & -93 & -63 & -1 & 0 & 0 & 0\\-293 & -868 & -1 & 0 & 0 & 0 & 879 & 2604 & 3\\0 & 0 & 0 & -293 & -868 & -1 & -1758 & -5208 & -6\\-1207 & -998 & -1 & 0 & 0 & 0 & 8449 & 6986 & 7\\0 & 0 & 0 & -1207 & -998 & -1 & -4828 & -3992 & -4\\-1218 & -309 & -1 & 0 & 0 & 0 & 3654 & 927 & 3\\0 & 0 & 0 & -1218 & -309 & -1 & 2436 & 618 & 2\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\end{matrix}\right]$$

Invert matrix $$P$$ and multiply with homogeneous vector to get matrix $$H$$

$$H = P ^{-1} * \displaystyle \left[\begin{matrix}0\\0\\0\\0\\0\\0\\0\\0\\1\end{matrix}\right]$$

Fill in and calculate result

$$H= \displaystyle \left[\begin{matrix}\frac{13724393133269}{1982067999646851}\\\frac{28367198390048}{1982067999646851}\\- \frac{5924445050578537}{660689333215617}\\\frac{11075363038922}{1982067999646851}\\- \frac{29759444826748}{1982067999646851}\\\frac{281612087155126}{660689333215617}\\\frac{2748231707}{1982067999646851}\\\frac{1890176860316}{1982067999646851}\\1\end{matrix}\right]$$

Convert to $$3 \times 3$$ matrix

$$H= \displaystyle \left[\begin{matrix}\frac{13724393133269}{1982067999646851} & \frac{28367198390048}{1982067999646851} & - \frac{5924445050578537}{660689333215617}\\\frac{11075363038922}{1982067999646851} & - \frac{29759444826748}{1982067999646851} & \frac{281612087155126}{660689333215617}\\\frac{2748231707}{1982067999646851} & \frac{1890176860316}{1982067999646851} & 1\end{matrix}\right]$$

In decimal numbers

$$H= \displaystyle \left[\begin{matrix}0.0069243 & 0.014312 & -8.9671\\0.0055878 & -0.015014 & 0.42624\\1.3865 \cdot 10^{-6} & 0.00095364 & 1.0\end{matrix}\right]$$

Check

General Formula to calculate point with the homography matrix

$$\left[\begin{array}{c}{x^{\prime} * \lambda} \\ {y^{\prime} * \lambda} \\ {\lambda}\end{array}\right]=\left[\begin{array}{lll}{h_{11}} & {h_{12}} & {h_{13}} \\ {h_{21}} & {h_{22}} & {h_{23}} \\ {h_{31}} & {h_{32}} & {h_{33}}\end{array}\right] \cdot\left[\begin{array}{l}{x} \\ {y} \\ {1}\end{array}\right]$$

Check with point $$p_5$$ $$(6|1|0)$$ in 3d, and $$(1540|502)$$ in 2d

$$\displaystyle \left[\begin{matrix}0.00692 & 0.0143 & -8.97\\0.00559 & -0.015 & 0.426\\1.39 \cdot 10^{-6} & 0.000954 & 1.0\end{matrix}\right] * \displaystyle \left[\begin{matrix}1540\\502\\1\end{matrix}\right] = \displaystyle \left[\begin{matrix}8.8809\\1.4942\\1.4809\end{matrix}\right]$$

$$x^{\prime} = 8.8809 / 1.4809 \approx 6$$

$$y^{\prime} = 1.4942 / 1.4809 \approx 1$$

Check with point $$p_6$$ $$(4|7|0)$$ in 3d, and $$(1679|128)$$ in 2d

$$\displaystyle \left[\begin{matrix}0.00692 & 0.0143 & -8.97\\0.00559 & -0.015 & 0.426\\1.39 \cdot 10^{-6} & 0.000954 & 1.0\end{matrix}\right] * \displaystyle \left[\begin{matrix}1679\\128\\1\end{matrix}\right] = \displaystyle \left[\begin{matrix}4.4907\\7.8863\\1.1244\end{matrix}\right]$$

$$x^{\prime} = 4.4907 / 1.1244 \approx 4$$

$$y^{\prime} = 7.8863 / 1.1244 \approx 7$$

SVD

Also described in the same linked question https://math.stackexchange.com/a/1289595/741822, the last Vector of $$V$$ of the SVD result should be the homography matrix. Unfortunately, when I calculate it, I get something which doesnt compare to the outcome before.

Im using a python script and the svd_r function described here http://mpmath.org/doc/current/matrices.html to do the job

from sympy import *
from mpmath import mp
from IPython.display import display

x_1 = [93,-7]
y_1 = [63,0]
x_2 = [293,3]
y_2 = [868,-6]
x_3 = [1207,7]
y_3 = [998,-4]
x_4 = [1218,3]
y_4 = [309,2]

P = Matrix([
[x_1[0]*-1 ,y_1[0]*-1 ,-1 ,0 ,0 ,0 ,x_1[0]*x_1[1] ,y_1[0]*x_1[1] ,x_1[1]],
[0 ,0 ,0 ,x_1[0]*-1 , y_1[0]*-1 ,-1 ,x_1[0]*y_1[1] ,y_1[0]*y_1[1] ,y_1[1]],
[x_2[0]*-1 ,y_2[0]*-1 ,-1 ,0 ,0 ,0 ,x_2[0]*x_2[1] ,y_2[0]*x_2[1] ,x_2[1]],
[0 ,0 ,0 ,x_2[0]*-1 , y_2[0]*-1 ,-1 ,x_2[0]*y_2[1] ,y_2[0]*y_2[1] ,y_2[1]],
[x_3[0]*-1 ,y_3[0]*-1 ,-1 ,0 ,0 ,0 ,x_3[0]*x_3[1] ,y_3[0]*x_3[1] ,x_3[1]],
[0 ,0 ,0 ,x_3[0]*-1 , y_3[0]*-1 ,-1 ,x_3[0]*y_3[1] ,y_3[0]*y_3[1] ,y_3[1]],
[x_4[0]*-1 ,y_4[0]*-1 ,-1 ,0 ,0 ,0 ,x_4[0]*x_4[1] ,y_4[0]*x_4[1] ,x_4[1]],
[0 ,0 ,0 ,x_4[0]*-1 , y_4[0]*-1 ,-1 ,x_4[0]*y_4[1] ,y_4[0]*y_4[1] ,y_4[1]],
])

U, S, V = mp.svd_r(P)
display(U)
display(S)
display(V)


which gives me

$$\displaystyle \left[\begin{matrix}0.0873968428008291 & 0.0682224411681068 & 7.78775669283991 \cdot 10^{-5} & -0.0307611970870902 & -0.047164344898438 & -4.31617782272004 \cdot 10^{-5} & -0.734415741368853 & -0.667210380053951 & -0.000763782447496494\\-0.1843309239413 & -0.0151678791434499 & -8.61313928411903 \cdot 10^{-5} & -0.177510439386387 & -0.196972727023475 & -0.000286352439937897 & 0.63453250533834 & -0.702034214967351 & -0.000453769669195275\\0.178502576804415 & 0.311156872128848 & 0.000325321190586577 & 0.77104600249753 & 0.446971881836657 & 0.000626740830458092 & 0.180594126282578 & -0.210730877512211 & -0.000435956847243289\\-0.669632242817859 & -0.535613720571974 & -0.000939878794635243 & 0.478237602951333 & -0.137967801791192 & 0.000272230101249859 & -0.129682866265852 & -0.0120304219769976 & -0.000973739592376422\\0.686925430229282 & -0.552023064176252 & 0.000248883448694341 & 0.254558652994286 & -0.38595799679894 & -0.000493796526964928 & 0.0865993870367576 & -0.046238819234477 & -0.00158707187572933\\-0.0790313367881196 & 0.55378718349052 & 0.0023403071836456 & 0.28199111075643 & -0.768889925477268 & -3.75757014301404 \cdot 10^{-5} & -0.0327206087187983 & 0.123639937181458 & 0.00038524047574376\\0.000450569782323263 & -0.0016191693558424 & 0.111042503882465 & 0.000924901054479592 & -0.000202764323620636 & 0.00454096268914545 & -0.00017046030825094 & -0.00107219827445236 & 0.993802818493512\\0.000319779779068856 & -0.000370212652047504 & 0.0463933646371045 & -0.000563323598922266 & -0.000441631139009571 & 0.998875291759086 & 0.000114547728791255 & -0.000112905332478902 & -0.0097483164699673\end{matrix}\right]$$

No data in that matrix compares to the calculated homography matrix above.

Creating a matrix with the last column and doing the check will end in $$\approx (0|0)$$

$$\displaystyle \left[\begin{matrix}-0.000763782447496494 & -0.000453769669195275 & -0.000435956847243289\\-0.000973739592376422 & -0.00158707187572933 & 0.00038524047574376\\0.993802818493512 & -0.0097483164699673 & 1\end{matrix}\right] * \displaystyle \left[\begin{matrix}1679\\128\\1\end{matrix}\right] = \displaystyle \left[\begin{matrix}-1.3409\\-1.8377\\1668.3\end{matrix}\right]$$

What am I doing wrong here?

Am I missing a step when working with SVD or is it that the 1st result (via the trick) just works by pure chance?

• I did not understand what you mean by 3D data and 2D data. A homography can only be calculated from points having similar dimensions. Jan 16, 2020 at 16:54
• Also, I think you need the last column of $V^T$, so take the last row of $V$. Jan 16, 2020 at 16:58
• I don't know what's going on with mpmath's svd function but when I try your code with numpy's svd it works just fine. Jan 16, 2020 at 17:27

As I pointed out in the comments, you need the last column of $$V^T$$, so take the last row of $$V$$.

Your code works just fine with numpy's svd as given below.

import numpy as np

x_1 = [93,-7]
y_1 = [63,0]
x_2 = [293,3]
y_2 = [868,-6]
x_3 = [1207,7]
y_3 = [998,-4]
x_4 = [1218,3]
y_4 = [309,2]

P = np.array([
[-x_1[0], -y_1[0], -1, 0, 0, 0, x_1[0]*x_1[1], y_1[0]*x_1[1], x_1[1]],
[0, 0, 0, -x_1[0], -y_1[0], -1, x_1[0]*y_1[1], y_1[0]*y_1[1], y_1[1]],
[-x_2[0], -y_2[0], -1, 0, 0, 0, x_2[0]*x_2[1], y_2[0]*x_2[1], x_2[1]],
[0, 0, 0, -x_2[0], -y_2[0], -1, x_2[0]*y_2[1], y_2[0]*y_2[1], y_2[1]],
[-x_3[0], -y_3[0], -1, 0, 0, 0, x_3[0]*x_3[1], y_3[0]*x_3[1], x_3[1]],
[0, 0, 0, -x_3[0], -y_3[0], -1, x_3[0]*y_3[1], y_3[0]*y_3[1], y_3[1]],
[-x_4[0], -y_4[0], -1, 0, 0, 0, x_4[0]*x_4[1], y_4[0]*x_4[1], x_4[1]],
[0, 0, 0, -x_4[0], -y_4[0], -1, x_4[0]*y_4[1], y_4[0]*y_4[1], y_4[1]],
])

[U, S, Vt] = np.linalg.svd(P)
homography = Vt[-1].reshape(3, 3)
transformedPoint = homography @ np.array([1679,  128, 1]).transpose()
print(transformedPoint/transformedPoint[-1]) # will be ~[4, 7, 1]
`
• this gives me a matrix $\displaystyle \left[\begin{matrix}0.000766575413886898 & 0.00158444869915194 & -0.992728903769269\\0.000618613946938351 & -0.00166221256660159 & 0.0471882946306038\\1.53502369083097 \cdot 10^{-7} & 0.000105575750777773 & 0.110708326585135\end{matrix}\right]$ which works, but is different tough. I thought they must be exactly the same. Jan 20, 2020 at 18:16
• Homographies are defined upto scale. This means given a homography matrix $H$, all homography matrices $kH$ where $k$ is a non-zero scalar are equivalent. In this case, if you multiply this matrix by $k = 9.02$ you will get the other matrix. Jan 20, 2020 at 22:25