initUndistortRectifyMap does not return same undistorted image as manual undistortion

Python
, , , ,
1

I have been trying to replicate what initUndistortRectifyMap does for my application. For some reason, when I do manual undistortion, I get different result than using initUndistortRectifyMap.

I feel I am following the steps to the letter, but somehow I get different outcome. Any insights will be greatly appreciated

Here is the undistorted image applying manually the undistortion equations (don’t mind the white lines in the picture, those are missing pixels that I didn’t interpolate for code simplicity). The forum doesn’t allow me to show two images in the same post but initUndistortRectifyMap-undistorted image looks more pincushiony for my camera parameters.

image
image3840×2160 362 KB

Here’s the code I used to repro the image above

import numpy as np
import cv2 as cv

# Generating an image of 2160x3840 pixels
# I am using the Charuco board generator that comes with OpenCV,
# (it can be any arbitrary image of 2160x3840 pixels)
num_squares_x = 32
num_squares_y = 18
square_width_m = 0.05 
marker_width_m = 0.04
dictionary = cv.aruco.getPredefinedDictionary(cv.aruco.DICT_6X6_1000)
board = cv.aruco.CharucoBoard(
   (num_squares_x, num_squares_y),
   square_width_m, marker_width_m, dictionary
)
board.setLegacyPattern(True)
img = board.generateImage((num_squares_x*120, num_squares_y*120))

h, w = img.shape[:2]

# Camera intrinsic parameters
f = 2931.87670606
cx = 1923.04908337
cy = 1070.12953259

# Distortion parameters (radial distortion only for my application)
k1 = 0.12027473
k2 = -0.49076343
k3 = 0.66055698
p1 = 0.00000000
p2 = 0.00000000

# my camera has same fx and fy
K = np.array([
    [  f,  0.0,   cx],
    [0.0,    f,   cy],
    [0.0,  0.0,  1.0]], dtype=np.float32)

dist = np.array([[k1, k2, p1, p2, k3]], dtype=np.float32)

# Method 1: using initUndistortRectifyMap (barrel)
map1, map2 = cv.initUndistortRectifyMap(K, dist, None, K, img.shape[:2][::-1], cv.CV_32FC1)
undistorted_imgv1 = cv.remap(img, map1, map2, cv.INTER_LINEAR)

# Method2 :  doing manual undistortion (pincushion)
#           (skiping pixel interpolation for simplicity)
undistorted_imgv2 = np.ones((h, w, 3), dtype=np.uint8) * 255

# Generating all the undistorted u',v' pixel coordinates from
# the original u,v pixel coordinates, using the camera intrinsics
# for the image with dimensions (h, w)
uv = np.meshgrid(range(0, w), range(0, h))
uu = uv[0].flatten()
vv = uv[1].flatten()

for i in range(len(uu)):
  u = uu[i]
  v = vv[i]
  xp = (u - cx) / f
  yp = (v - cy) / f
  r2 = xp**2 + yp**2
  r4 = r2**2
  r6 = r2**3
  gamma = 1 + k1*r2+ k2*r4+ k3*r6
  xpp = xp * gamma + 2*p1*xp*yp + p2*(r2 + 2*xp**2)
  ypp = yp * gamma + p1*(r2+ 2*yp**2) + 2*p2*xp*yp
  new_uv=np.dot(K, np.array([xpp, ypp, 1]).reshape(3,1)).flatten()
  up = new_uv[0] / new_uv[2]
  vp = new_uv[1] / new_uv[2]
  # Check if the new pixel coordinates are within the bounds of the image
  if up < 0 or up >= w or vp < 0 or vp >= h:
    continue
  undistorted_imgv2[int(vp), int(up)] = img[int(v), int(u)]

cv.imwrite("original_image.png", img)
cv.imwrite("undistorted_image_initUndistortRectifyMap.png", undistorted_imgv1)
cv.imwrite("undistorted_manual.png", undistorted_imgv2)
2

Here’s how initundistortrectifymap-undistorted image looks like (more pincushion)

image
image3840×2160 216 KB

3

this the original image generated with the code above

image
image3840×2160 65.2 KB

4

looks to me like

  1. appearance: the maps are opposites/inverses of each other (picture corners dragged inward/outward)
  2. image mapping: you pushed pixels from a source to a destination space (assignment into undistorted_imgv2)

some basic statements:

you could, theoretically, model the lens with either a distortion model or an undistortion model.

these models are complex enough (polynomials at least) to be impossible/intractable to invert algebraically. they require numerical inversion. that’s a few iterations of simple arithmetic at least, so a couple times more expensive than a single evaluation.

numerical inversion of a few points is cheap enough. numerical inversion of an entire image gets costly. hence the choice of equations in such a direction to make image mapping cheap (i.e. doesn’t require inversion, only the direct evaluation of the model equation).

mapping of images: you pushed (2) when you should have pulled. the pushing causes those moire patterns (lines), where density varies due to the map. yes, they’re moire patterns. it’s the same principle. when pulling, as is done by cv::remap() and all the warp functions in OpenCV, you consider each destination pixel (hence uniform density), map its position into the source space, and sample the source at that location. that calculation requires an inverse map to what you would consider the “forward” direction. mapping an image requires a “backwards” equation. if you called one of the warp functions (affine, perspective), you would give it a single forward transform matrix to represent that map, and the function internally inverts it (or numerically sensible equivalent operation).

in practice, people want undistortion of images to be fast (and undistortion of points to be cheap enough). that motivates the choice of an “inverse undistortion” model, which cancels to a distortion model. the model being forward (distortion) has the added benefit of being nice to think about. you can think about ideal rays or points in space, and what the lens does to them, instead of image points and where they came from due to the lens distortion.

if you need to generate distorted views, there is a reasonably simple and cheap way to generate a distorting (inverse undistorting) map for cv::remap() by way of numerically inverting the undistorting (inverse distorting) map gained from initUndistortRectifyMap. https://stackoverflow.com/questions/41703210/inverting-a-real-valued-index-grid

I don’t know if OpenCV contains a specific function for such a thing. it might. I never specifically needed to calculate distorted images, only undistorted images or distorted points.