Homography transformations describe the change in projected position of points on a rigidly moving plane observed between 2D images.
A homography H ∈ ℝ³ˣ³ transforms a point (x y)∈ℝ² to a new point (x' y')∈ℝ² by linear transformation in homogeneous coordinates:
/ x'/w \ / h₁₁ h₁₂ h₁₃ \ / xᵢ \
| y'/w | = | h₂₁ h₂₂ h₂₃ | | yᵢ |
\ w / \ h₃₁ h₃₂ h₃₃ / \ 1 /
We can pose the problem of finding the homography H which matches a set of input points (xᵢ yᵢ)∈ℝ² to corresponding target points (x'ᵢ y'ᵢ)∈ℝ²:
‖ / x'ᵢ \ - / (h₃₁ h₃₂ h₃₃ )/ xᵢ \ \⁻¹ /h₁₁ h₁₂ h₁₃ \/ xᵢ \ ‖²
min ∑ᵢ ‖ \ y'ᵢ / | | yᵢ | | \h₂₁ h₂₂ h₂₃ /| yᵢ | ‖
H ‖ \ \ 1 / / \ 1 / ‖
The term inside the square is non-linear in H. The common adjustment is to multiply this term through so that the result is a linear least-squares problem:
‖ / (h₃₁ h₃₂ h₃₃ )/ xᵢ \ \ / x'ᵢ \ - /h₁₁ h₁₂ h₁₃ \/ xᵢ \ ‖²
min ∑ᵢ ‖ | | yᵢ | | \ y'ᵢ / \h₂₁ h₂₂ h₂₃ /| yᵢ | ‖
H ‖ \ \ 1 / / \ 1 / ‖
This can be rearranged as:
‖ / -xᵢ -yᵢ -1 0 0 0 x'ᵢxᵢ - x'ᵢyᵢ -x'ᵢ \ / h₁₁ \ ‖²
‖ \ 0 0 0 -xᵢ -yᵢ -1 y'ᵢxᵢ - y'ᵢyᵢ -y'ᵢ / | h₁₂ | ‖
‖ ^---------------------.---------------------^ | h₁₃ | ‖
‖ Aᵢ | h₂₁ | ‖
min ∑ᵢ ‖ | h₂₂ | ‖
H ‖ | h₂₃ | ‖
‖ | h₃₁ | ‖
‖ | h₃₂ | ‖
‖ \ h₃₃ / ‖
and further simplified as:
min (h₁₁ h₁₂ h₁₃ h₂₁ h₂₂ h₂₃ h₃₁ h₃₂ h₃₃) (∑ᵢ AᵢᵀAᵢ) / h₁₁ \
H ^----Q---^ | h₁₂ |
| h₁₃ |
| h₂₁ |
| h₂₂ |
| h₂₃ |
| h₃₁ |
| h₃₂ |
\ h₃₃ /
To avoid the trivial solution H=0, we take the first eigen vector of Q, which corresponds to the minimizer subject to the constraint H:H = 1. Other sources suggest constraining h₃₃=1.
Unfortunately, this solution does not minimize the original least-squares problem. The trouble started when we multiplied the term inside the square. This can be understood instead by introducing a per-point bias wᵢ²:
‖ / x'ᵢ \ - / (h₃₁ h₃₂ h₃₃ )/ xᵢ \ \⁻¹ /h₁₁ h₁₂ h₁₃ \/ xᵢ \ ‖²
min ∑ᵢ wᵢ² ‖ \ y'ᵢ / | | yᵢ | | \h₂₁ h₂₂ h₂₃ /| yᵢ | ‖
‖ \ \ 1 / / \ 1 / ‖
where
wᵢ = (h₃₁ h₃₂ h₃₃ )/ xᵢ \
| yᵢ |
\ 1 /
Perhaps without noticing it, we've placed more weight on reproducing points with a large homogeneous coordinate. The more perspective distortion in the observation, the stronger this bias will be.
Alternatively, we could simply minimize the original least-squares problem directly using Gauss-Newton method:
‖ / x'ᵢ \ - / (h₃₁ h₃₂ h₃₃ )/ xᵢ \ \⁻¹ /h₁₁ h₁₂ h₁₃ \/ xᵢ \ ‖²
min ∑ᵢ ‖ \ y'ᵢ / | | yᵢ | | \h₂₁ h₂₂ h₂₃ /| yᵢ | ‖
H ‖ \ \ 1 / / \ 1 / ‖
^----------------------------Fᵢ(H)----------------------^
Given some current guess H we can define the per-point Jacobians Jᵢ as:
Jᵢ = ∂Fᵢ/∂H = (tᵢ zᵢ) ∈ ℝ²ˣ⁹
tᵢ = / (h₃₁ h₃₂ h₃₃ )/ xᵢ \ \⁻¹ / xᵢ yᵢ 1 0 0 0 \
| | yᵢ | | \ 0 0 0 xᵢ yᵢ 1 /
\ \ 1 / /
zᵢ = / (h₃₁ h₃₂ h₃₃ )/ xᵢ \ \⁻² / xᵢ yᵢ 1 \ ⊙ /h₁₁ h₁₂ h₁₃ \/ xᵢ \
| | yᵢ | | \ xᵢ yᵢ 1 / \h₂₁ h₂₂ h₂₃ /| yᵢ |
\ \ 1 / / \ 1 /
Then the Gauss-Newton update is written as:
H ← H - (∑ᵢ JᵢᵀJᵢ)⁺ (∑ᵢ Jᵢᵀ Fᵢ)
where M⁺ is the pseudo-inverse of M.
It seems sufficient to initialize H = identity and take just 3 or 4 Gauss-Newton steps, implying that this should just be a small constant factor slower than non-optimal Eigen-decomposition method above. In many problems, when a better initial guess may be available, a single step may be enough.
It's not immediately clear to me whether the original least squares problem is convex over H. If not, then there should be some concern of finding a local rather than global minimum.