Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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 vectors starts $(x,y)$ and ends $(x',y')$ and I know the distortion equations looks like this: $$ 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 determined through measurements, they are most likely biased.

How can I calculate $k_1$ and $k_2$?

share|improve this question

2 Answers 2

up vote 2 down vote accepted

Pick two starting points $(x_1,y_1)$ and $(x_2,y_2)$ that are not colinear 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^b$.

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

share|improve this answer
    
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? –  streuguut Feb 13 '13 at 17:36
    
Ok, makes sense. I've explained an approach for that in the last edit. –  Elmar Zander Feb 13 '13 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. –  Elmar Zander Feb 13 '13 at 19:12
    
OTOH: I think it's more readable now, somehow. –  Elmar Zander Feb 13 '13 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! –  streuguut Feb 14 '13 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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.