Calculate Homography with and without SVD

Question

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 doesn`t 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. — ImmenselyHappy, Jan 16 '20 at 16:54
Also, I think you need the last column of $V^T$, so take the last row of $V$. — ImmenselyHappy, Jan 16 '20 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. — ImmenselyHappy, Jan 16 '20 at 17:27

score 5 · Accepted Answer · edited Dec 07 '20 at 22:07

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. — hypnomaki, Jan 20 '20 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. — ImmenselyHappy, Jan 20 '20 at 22:25

Calculate Homography with and without SVD

3D Data

2D Data

Calculation

Check

SVD

1 Answers1