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?
The opencv undistortion may be a bit similar, so that link may be useful but it is not quite the same though.
