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

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

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.

• OK, I see. But to take my question one step further. If I got my distorted points only through measurements, then I should use all of them and put them into some sort of an optimization equation, wouldn't I?
– Fei
Feb 13, 2013 at 17:36
• Ok, makes sense. I've explained an approach for that in the last edit. Feb 13, 2013 at 19:07
• @streuguut It's a bit unfortunate, that you changed your naming of variables. I'll try to reflect this in my answer. Feb 13, 2013 at 19:12
• OTOH: I think it's more readable now, somehow. Feb 13, 2013 at 19:18
• Yes, changing the notation made it easier to write down the QR decomp. part and I was already on my way to edit your answer, but my edits didn't pass the peer review. ^^" But now both the question and the answer are matching perfectly. Thanks!
– Fei
Feb 14, 2013 at 9:57

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.