I have this vector field full of displacement vectors, which indicates radial distortions by a lens system.
(Source)
I know where each of the displacement vectors starts $(x,y)$ and ends $(x',y')$ and I know the distortion equations look like
$$
x' = (1 + k_1r^2 + k_2r^4)x\\
y' = (1 + k_1r^2 + k_2r^4)y
$$
where $r^2 = x^2 + y^2$. Since each $(x',y')$ was measured, they are most likely biased. How can I estimate coefficients $k_1$ and $k_2$?
Best Answer
Pick two starting points $(x_1,y_1)$ and $(x_2,y_2)$ that are not collinear with the origin, and the corresponding end points. Then you get the two equations $$ x'_1 = (1+k_1r_1^2+k_2r_1^4)x_1\\ x'_2 = (1+k_1r_2^2+k_2r_2^4)x_2. $$ Then rewrite as a linear equation $$ \begin{bmatrix}x'_1/x_1-1\\x'_2/x_2-1\end{bmatrix}= \begin{bmatrix}r_1^2 & r_1^4 \\ r_2^2 & r_2^4 \end{bmatrix} \begin{bmatrix}k_1\\k_2\end{bmatrix} $$ and solve for $k_1$ and $k_2$.
Edit: If you want to include all measured vectors into your estimation of $k_1$ and $k_2$, you can set up a linear system like $Ax=b$ follows: $$Ax= \begin{bmatrix}r_1^2 & r_1^4 \\r_1^2 & r_1^4 \\ r_2^2 & r_2^4 \\ r_2^2 & r_2^4 \\ \vdots & \vdots \\ r_n^2 & r_n^4 \end{bmatrix} \begin{bmatrix}k_1\\k_2\end{bmatrix}= \begin{bmatrix}x'_{1}/x_{1}-1\\y'_{1}/y_{1}-1\\x'_{2}/x_{2}-1\\y'_{2}/y_{2}-1\\\vdots\\y'_{n}/y_{n}-1\end{bmatrix} =b$$ This system is clearly overdetermined ($A$ is a $2n\times2$ matrix), so you need to employ least squares methods to solve it. Two methods to solve it:
Normal equations: premultiply both sides of the equation by $A^T$ and solve $A^TAx=A^Tb$.
QR decomposition: find a decomposition $QR=A$ such that $Q\in\mathbb R^{2\times n}$ has orthogonal columns, i.e. $Q^TQ=I$, and $R$ is an upper triangular $2\times2$ matrix. Then solve $Rx=Q^Tb$ (to see this just premultiply $Ax=b$ by $Q^T$).
The second method is numerically more stable, and you'll find the $QR$ decomposition in nearly all linear algebra packages or mathematical software environments. If you have Matlab, you can directly use
b=A\x, since matlab will automatically use a least squares method if $A$ is not square.