Approximating the inverse of a polynomial function for camera calibration

Question

The method below to remove lens distortion from a camera was written more than ten years ago and i am trying to understand how the approximation works.

void Distortion::removeLensDistortion(const V2 &dist, V2 &ideal) const
{
 const int ITERATIONS=10;
 V2 xn;
 xn.x=(dist.x - lensParam.centerX)/lensParam.fX;
 xn.y=(dist.y - lensParam.centerY)/lensParam.fY;

 V2 x=xn;
 for (int i=0;i<ITERATIONS;i++) {
   double r2 = Utilities::square(x.x)+Utilities::square(x.y);
   double r4 = Utilities::square(r2);
   double rad=1+lensParam.kc1 * r2 + lensParam.kc2 * r4;
   x.x/=rad;
   x.y/=rad;
 }
 ideal.x=x.x*lensParam.fX+lensParam.centerX;
 ideal.y=x.y*lensParam.fY+lensParam.centerY;
}

As a reminder:

• lensParam.centerX and lensParam.centerY is the principal point
• lensParam.fX and lensParam.fY is the focal length in pixel
• lensParam.kc1 and lensParam.kc2 are the first two radial distortion coefficients. This is $k_1$ and $k_2$ in the formula below.

The formula to add lens distortion given the first two radial distortion parameters is as follows:

$$ x_{distorted} = x_{undistorted} \cdot (1+k_1\cdot r^2 + k_2 \cdot r^4) \\ y_{distorted} = y_{undistorted} \cdot (1+k_1\cdot r^2 + k_2 \cdot r^4) $$ where $r^2=x_{undistorted}^2+y_{undistorted}^2$ and $r^4=(r^2)^2$

In the code above, the term $(1+k_1\cdot r^2 + k_2 \cdot r^4)$ is calculated and saved in the variable rad and the distorted x is divided by rad in each of the ten iterations. All of the cameras we use have a pincushion distortion (so $k_1<0$)

The question is how is this algorithm approximating the undistorted image points? Do you know if there is any paper in which this algorithm is proposed?

This is not a general question about approximating the inverse of polynomials but rather in the context of camera calibration. However, it may still boil down to pure math. Tag suggestions would be helpful.

Thank you in advance.

The opencv undistortion may be a bit similar, so that link may be useful but it is not quite the same though.

Answer 1

Warning: my answer is at odds with your code (my C++ is rusty)

You are given

$x_{d}=x_{u}F(x_u,y_u)$

$y_{d}=y_{u}F(x_u,y_u)$

where the $\cdot _d,\cdot_u$ refer to distorted/undistorted respectively. However you want to solve for $x_u;y_u$ given $x_d,y_d$. Write

$x_{u}={x_{d}\over F(x_u,y_u)}$

$y_{u}={y_{d}\over F(x_u,y_u)}$

but you don't know the values of $x_u;y_u$ are, so you solve by iteration:

$x^{n+1}_{u}={x_{d}\over F(x^n_u,y^n_u)}$

$y^{n+1}_{u}={y_{d}\over F(x^n_u,y^n_u)}$

and start with an initial guess for $x_u$ to be $x^0_u = x_d$ (and similarly for $y_u$ )

However that is NOT what the code seems to be doing. Your code implements the iteration $x^{n+1}_{u}={x^n_{u}\over F(x^n_u,y^n_u)}$. So something is wrong in my explanation or your code. I am leaving the answer here, and if it turns out I am wrong I will delete it.