# Numerically Stable Homography Computation

Suppose you want to find a Homography transform that maps a set of 2D coordinates to another arbitrary set of 2D coordinates. This problem has many applications, such as image registration based on keypoint matching. This post will describe a few steps to make sure the homography estimation is numerically stable and can be performed using 32bit floating point coordinates.

## Mathematical Solution

Firstly, we should formalize the problem definition. Let $C=\{\{p_1,p^\prime_1\},\{p_2,p^\prime_2\}, \dots,\{p_n,p^\prime_n\}\}, \; n \geq 4$, be a set of n pairs of 2D coordinates for which there exists a homography transform $H$ such that $H(p_i)=p^\prime_i$ for all $\{p_i, p^\prime_i\} \in C$. The transform equation can be written in matrix form as below:

$\langle p^{\prime}_{xi}, p^{\prime}_{yi} \rangle=\langle \frac{h_1 p_{xi} + h_2 p_{yi} + h_3}{h_7 p_{xi} + h_8 p_{yi} + 1}, \frac{h_4 p_{xi} + h_5 p_{yi} + h_6}{h_7 p_{xi} + h_8 p_{yi} + 1} \rangle,$

where $h_1, \dots, h_8$ are the coefficients of the homography transform. This equation can be rewritten in matrix form as

$\begin{bmatrix} p_{x1}& p_{y1}& 1 & 0 & 0 & 0 & -p^\prime_{x1}p_{x1} & -p^\prime_{x1}p_{y1} & -p^\prime_{x1} \\ 0 & 0 & 0 & p_{x1}& p_{y1}& 1 &-p^\prime_{y1}p_{x1} & -p^\prime_{y1}p_{y1} & -p^\prime_{y1} \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ p_{xn}& p_{yn}& 1 & 0 & 0 & 0 & -p^\prime_{xn}p_{xn} & -p^\prime_{xn}p_{yn} & -p^\prime_{xn} \\ 0 & 0 & 0 & p_{xn}& p_{yn}& 1 &-p^\prime_{yn}p_{xn} & -p^\prime_{yn}p_{yn} & -p^\prime_{yn} \end{bmatrix} \begin{bmatrix} h_1 \\ h_2 \\ h_3 \\ h_4 \\ h_5 \\ h_6 \\ h_7 \\ h_8 \\ 1 \end{bmatrix} =\mathbf{0}. \; (I)$

This is a linear system in the form $Ah=\mathbf{0}$, with $A$ having 2n rows and 9 columns. To avoid the pitfall of a trivial solution (which would render $h=\mathbf{0}$), this system can be solved by using the PCA algorithm, as follows.

If we multiply both sides of equation (I) by $A^t$, we end up with an equation in the form $A^tAh=\mathbf{0}$. Note that $A^tA$ is a square matrix of size 8×8. Therefore, it is possible to compute its eigenvalues and eigenvectors. The final value of $h$ is the eigenvector of $A^tA$ with the lowest eigenvalue. Note that we may need to normalize the final solution by dividing all elements of $h$ by its eighth element.

## Numerical Solution

As we have seen, in order to determine the homography parameters, we need to calculate the matrix $A^tA$. It is easy to verify that the largest element of $A^tA$ may scale up to quartically with respect to the input coordinates, while its smallest elements scale linearly with respect to the input coordinates.

$(A^tA)(1:6, :)=\\\begin{bmatrix} \sum\limits^n_{i=1} p^2_{xi} & \sum\limits^n_{i=1} p_{xi}p_{yi} & \sum\limits^n_{i=1} p_{xi} & 0 & 0 & 0 & \sum\limits^n_{i=1}-p^\prime_{xi}p^2_{xi} & \sum\limits^n_{i=1}-p^\prime_{xi}p_{xi}p_{yi} & \sum\limits^n_{i=1}-p^\prime_{xi}p_{xi} \\ \sum\limits^n_{i=1} p_{xi}p_{yi} & \sum\limits^n_{i=1} p^2_{yi} & \sum\limits^n_{i=1} p_{yi} & 0 & 0 & 0 & \sum\limits^n_{i=1}-p^\prime_{xi}p_{xi}p_{yi} & \sum\limits^n_{i=1}-p^\prime_{xi}p^2_{yi} & \sum\limits^n_{i=1}-p^\prime_{xi}p_{yi} \\ \sum\limits^n_{i=1} p_{xi} & \sum\limits^n_{i=1} p_{yi} & n & 0 & 0 & 0 & \sum\limits^n_{i=1}-p^\prime_{xi}p_{xi} & \sum\limits^n_{i=1}-p^\prime_{xi}p_{yi} & \sum\limits^n_{i=1}-p^\prime_{xi} \\ 0 & 0 & 0 & \sum\limits^n_{i=1} p^2_{xi} & \sum\limits^n_{i=1} p_{xi}p_{yi} & \sum\limits^n_{i=1} p_{xi} & \sum\limits^n_{i=1}-p^\prime_{yi}p^2_{xi} & \sum\limits^n_{i=1}-p^\prime_{yi}p_{xi}p_{yi} & \sum\limits^n_{i=1}-p^\prime_{yi}p_{xi} \\ 0 & 0 & 0 & \sum\limits^n_{i=1} p_{xi}p_{yi} & \sum\limits^n_{i=1} p^2_{yi} & \sum\limits^n_{i=1} p_{yi} & \sum\limits^n_{i=1}-p^\prime_{yi}p_{xi}p_{yi} & \sum\limits^n_{i=1}-p^\prime_{yi}p^2_{yi} & \sum\limits^n_{i=1}-p^\prime_{yi}p_{yi} \\ 0 & 0 & 0 & \sum\limits^n_{i=1} p_{xi} & \sum\limits^n_{i=1} p_{yi} & n & \sum\limits^n_{i=1}-p^\prime_{yi}p_{xi} & \sum\limits^n_{i=1}-p^\prime_{yi}p_{yi} & \sum\limits^n_{i=1}-p^\prime_{yi} \end{bmatrix}$

As an example, consider the first six rows of $A^tA$ above and compare them with the seventh row below, which was transposed for clarity:

$(A^tA)(7, :)^t = \\ \begin{bmatrix} \sum\limits^n_{i=1}-p^\prime_{xi}p^2_{xi} \\ \sum\limits^n_{i=1}-p^\prime_{xi}p_{xi}p_{yi} \\ \sum\limits^n_{i=1}-p^\prime_{xi}p_{xi} \\ \sum\limits^n_{i=1}-p^\prime_{yi}p^2_{xi} \\ \sum\limits^n_{i=1}-p^\prime_{yi}p_{xi}p_{yi} \\ \sum\limits^n_{i=1}-p^\prime_{yi}p_{xi} \\ \sum\limits^n_{i=1}(p^\prime_{yi}p_{xi})^2 \\ \sum\limits^n_{i=1}p^{\prime2}_{xi}p_{xi}p_{yi} \\ \sum\limits^n_{i=1}p^{\prime2}_{xi}p_{xi} \\ \end{bmatrix}.$

Thus, it is clear that we may end up with a matrix with wildly different numbers in terms of scale if the input coordinates are not preprocessed. This difference can be large enough for significant numerical imprecisions to arise if we mix them in the same computation, depending on the underlying primitive type used to represent floating point numbers. This, however, is exactly what happens in numerical methods for determining eigenvectors and eigenvalues.

Therefore, we need to cope with this issue if we want to obtain a precise homography estimate. One way to do so would be to use a broader primitive type to represent our matrix (e.g. use double instead of float), but that would only be mitigating the real problem. A better solution is to normalize the input coordinates so that all of them fit in the unit range $-1 \leq p,p^\prime \leq 1$.

This can be done by selecting a scale parameter $s$ such that $s=\max{p_{xi}, p_{yi}, p^\prime_{xi}, p^\prime_{yi}} \; \forall \; 1 \leq i \leq n$. Then, we divide all coordinates $p$ and $p^\prime$ by $s$ and compute the matrix $\tilde{A}$ using these normalized coordinates. Finally, we find the homography parameters for $\tilde{A}^t\tilde{A}$.

Note that this homography maps normalized coordinates to normalized coordinates, i.e. $\tilde{H}(p_i/s)=p^\prime_i/s$. To find the final homography that maps the unscaled coordinates, we need to perform the following composition:

$H=S^{-1} \circ \tilde{H} \circ S,$

where $S$ corresponds to the function that divides a coordinate by $s$.

## Final thoughts

To illustrate the idea of the discussed method for computing the homography parameters, consider the Matlab/Octave implementation below.

% Input coordinates in homogeneous coordinates
p=[
281.1662 516.9434 484.2327 262.9684
154.7470 136.7685 379.9645 379.7526
1        1        1        1
];
d=[
290 490 490 290
159 159 359 359
1   1   1   1
];

% Scale parameter
s=max(max(p(:), d(:)));

% Number of pairs of coordinates
number_coords = size(p, 2);

% Homography transform function
H=@(p, h) h*p ./ repmat((h*p)(3,:), 3, 1);

% Normalize coordinates
p_s = p / s;
d_s = d / s;

% Build matrix A
A=[];

for i=1:number_coords
A = [
A;
p_s(1,i) p_s(2,i) 1 ...
0      0      0 ...
-p_s(1,i)*d_s(1,i) -p_s(2,i)*d_s(1,i) -d_s(1,i);
0      0      0 ...
p_s(1,i) p_s(2,i) 1 ...
-p_s(1,i)*d_s(2,i) -p_s(2,i)*d_s(2,i) -d_s(2,i)
];
end

% Compute eigenvectors and eigenvalues
[v, l] = eig(A'*A);

% Normalize and convert homography to matrix form
h_s = reshape(v(:,1)/v(9,1), [3 3])';

% Represent scale as a matrix
S = [
1/s 0 0
0 1/s 0
0 0 1
];
S_inv = [
s 0 0
0 s 0
0 0 1
];

% Compose normalized homography with scale operators
h = S_inv * h_s * S;

% Test: d ~= H(p,h)
abs(H(p,h) - d)