How to calculate the lens distortion coefficients with a known displacement vector field?

I have this vector field full of displacement vectors, which indicates radial distortions by a lens system.

vector field

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

Solutions Collecting From Web of "How to calculate the lens distortion coefficients with a known displacement vector field?"

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:

  1. Normal equations: premultiply both sides of the equation by $A^T$ and solve $A^TAx=A^Tb$.

  2. 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.

For a single point $(x_1,y_1)$ and its conjugate $(x_2,y_2)$:

$$x_2^2 + y_2^2 = (1+k_1 r^2 + k_2 r^4)^2 (x_1^2 + y_1^2) = (1+k_1 r^2 + k_2 r^4)^2 r^2$$

so that

$$1+k_1 r^2 + k_2 r^4 = \pm \sqrt{\frac{x_2^2+y_2^2}{r^2}} $$

The sign of course depends on whether the distortion also inverts. In your case of simple pincushion distortion, I don’t think so, so use the positive solution.

Repeat this for another point and its conjugate and you have two equations and two unknowns which will determine $k_1$ and $k_2$. You may wish to do this over the pupil/field to get a sense of how well your imaging obeys this model.