Introduction

When someone is trying to get the internal parameter of a camera, this person thinks about calibrating the camera using a chessboard. Then one would naturally follow this tutorial, getting the results from the good opencv project. Everything is great.

After such calibration, you can undistort an image or features inside the image with the obtained distortion coefficients. But, there are two issues, (at least for me)

  1. What does the numbers in the distortion coefficients mean?
  2. How do I use these numbers to distort & undistort features without calling functions in opencv?

Not the Truth

Equations from opencv

Going back to the tutorial, one see these equations for correcting the ,

\[\begin{aligned} x_\text{corrected} &= x (1 + k_1 r^2 + k_2 r^4 + k_3 r^6) & \text{Radial Distortion}\\ y_\text{corrected} &= y (1 + k_1 r^2 + k_2 r^4 + k_3 r^6) \\[0.5em] x_\text{corrected} &= x + [2 p_1 x y + p_2(r^2 + 2x^2)] & \text{Tangential Distortion} \\ y_\text{corrected} &= y + [p_1(r^2 + 2y^2) + 2p_2 xy] \end{aligned}\]

However, when I implement above equations to undistort points, I get strange results. The following part is my implementation,

def undistort_wrong(xy, k, distortion):
    k1, k2, p1, p2, k3 = distortion
    fx, fy = k[0, 0], k[1, 1]
    cx, cy = k[:2, 2]
    x, y = xy.astype(float)
    x0 = (x - cx) / fx
    y0 = (y - cy) / fy
    r2 = x0 ** 2 + y0 ** 2
    k = (1 + k1 * r2 + k2 * r2**2 + k3 * r2**3)
    delta_x = 2 * p1 * x0*y0 + p2 * (r2 + 2 * x0**2)
    delta_y = p1 * (r2 + 2 * y0**2) + 2 * p2 * x0*y0
    x = x0 * k + delta_x
    y = y0 * k + delta_y
    return np.array((x * fx + cx, y * fy + cy))

And the result is what I get, and it is wrong. The tomato crosses are the undistorted points from my code and the teal circles are results from cv2.undistortPoints. They should match.

Matlab version

Another version is from the matlab tutorial, which is contrary to the opencv version, but looks correct

\[\begin{aligned} x_\text{distorted} &= x (1 + k_1 r^2 + k_2 r^4 + k_3 r^6) & \text{Radial Distortion}\\ y_\text{distorted} &= y (1 + k_1 r^2 + k_2 r^4 + k_3 r^6) \\[0.5em] x_\text{distorted} &= x + [2 p_1 x y + p_2(r^2 + 2x^2)] & \text{Tangential Distortion} \\ y_\text{distorted} &= y + [p_1(r^2 + 2y^2) + 2p_2 xy] \end{aligned}\]

It is confusing but these turned out to be the actual equations that opencv implemented.

The Truth

Source code in opencv

The truth is, as usual, in the source code of opencv. Typically, we see the relevant code blocks,

x0 = x = (x - cx)*ifx;
y0 = y = (y - cy)*ify;

// compensate distortion iteratively
for( j = 0; j < iters; j++ )
{
    double r2 = x*x + y*y;
    double icdist = (1 + ((k[7]*r2 + k[6])*r2 + k[5])*r2)/(1 + ((k[4]*r2 + k[1])*r2 + k[0])*r2);
    double deltaX = 2*k[2]*x*y + k[3]*(r2 + 2*x*x);
    double deltaY = k[2]*(r2 + 2*y*y) + 2*k[3]*x*y;
    x = (x0 - deltaX)*icdist;
    y = (y0 - deltaY)*icdist;
}

This confirmed my guess and the math equations in the opencv tutorial is not correct.

Undistort manually

The following code correctly undistorts points,

def undistort(xy, k, distortion, iter_num=3):
    k1, k2, p1, p2, k3 = distortion
    fx, fy = k[0, 0], k[1, 1]
    cx, cy = k[:2, 2]
    x, y = xy.astype(float)
    x = (x - cx) / fx
    x0 = x
    y = (y - cy) / fy
    y0 = y
    for _ in range(iter_num):
        r2 = x ** 2 + y ** 2
        k_inv = 1 / (1 + k1 * r2 + k2 * r2**2 + k3 * r2**3)
        delta_x = 2 * p1 * x*y + p2 * (r2 + 2 * x**2)
        delta_y = p1 * (r2 + 2 * y**2) + 2 * p2 * x*y
        x = (x0 - delta_x) * k_inv
        y = (y0 - delta_y) * k_inv
    return np.array((x * fx + cx, y * fy + cy))

And this is the result, it is correct!