Where are the inverse equations hiding?

A common problem in computer vision is modelling lens distortion.  Lens distortion distort the shape of objects and introduce large errors in structure from motion applications.  Techniques and tools for calibrating cameras and removing lens distortion are widely available.  While books and papers readily provide the forward distortion equations (given an ideal undistorted coordinate, compute the distorted coordinate) inverse equations are much harder to come by.

Turns out there is no analytic equation for inverse radial distortion.  Which might explain the mysterious absence of inverse equations, but still it would be nice if this issue was discussed in common references. First a brief summary of background equations is given, followed by the solution to the inverse distortion problem.

Inverse Distortion Problem: Given a distorted image coordinate and distortion parameters, determined the coordinate of the ideal undistorted coordinate.

Background

Unlike the idealized pin-hole camera model, real world cameras have lens distortions.  Radial distortion is a common model for lens distortion and is summarized by the following equations: $\hat{x} = x + x[k_1 r^2 + k_2 r^4]$ $\hat{u} = u + (u-u_0)[k_1 r^2 + k_2 r^4]$ $r^2=u_x^2 + u_y^2$

where $\hat{u}=[\hat{u}_x,\hat{u}_y]$ and $u$ are the observed distorted and ideal undistorted pixel coordinates, $\hat{x}$ and $x$ are the observed and ideal undistorted normalized pixel coordinates, and $k_1, k_2$ are radial distortion coefficients.  The relationship between normalized pixel and unnormalized pixel coordinates is determined by the camera calibration matrix, $u=C x$, where $C$ is the 3×3 calibration matrix.

Solution

While there is no analytic solution there is an iterative solution.  The iterative solution works by first estimating the radial distortion’s magnitude at the distorted point and then refining the estimate until it converges.

1. $x=C^{-1}\hat{u}$
2. do until converge:
1. $r=||x||^2$
2. $m=k_1 r^2 + k_2 r^4$
3. $x=\frac{C^{-1}\hat{u}}{1+m}$
3. $u = \frac{x+x_0 m}{1+m}$

Where $x_0$ is the principle point/image center, and $||x||$ is the Euclidean norm. Only a few iterations are required before convergence.  For added speed when applying to an entire image the results can be cached.  A similar solution can be found by adding in terms for tangential distortion.