Issue with Radial Distortion in Image Processing - Getting Barrel Distortion Instead of Pincushion Distortion

I'm encountering an issue with radial distortion in image processing. I'm trying to apply radial distortion to an image using distortion coefficients k_1 = 1 , k_2 = 0.4 and k_3 = 0.2 . According to the distortion model, I should be getting pincushion distortion as the output. However, when I tried two different methods, I got two different results.

Method 1: Bilinear Interpolation with Array Calculation

In this method, I used pure array calculation to apply radial distortion, the output is pincushion distortion. Here's the code I used:

def apply_distortion_jit(img,W,img_mat, k1, k2, k3):
    H = W
    dis_center_x = W/2
    dis_center_y = H/2

    for x in range(W-2):
        for y in range(H-2):
            x_norm = (x - dis_center_x)/ (W)
            y_norm = (y - dis_center_y)/ (H)
            r = x_norm * x_norm + y_norm * y_norm

            coeff = (k1 + k2 * r + k3 * (r ** 2))
            x_d = coeff * (x - W/2) + W/2
            y_d = coeff * (y - H/2) + H/2
            if x_d >= W-1 or y_d >= H-1 or x_d<0 or y_d <0:
                continue

            x_d_int = math.floor(x_d)
            y_d_int = math.floor(y_d)

            dx = x_d - x_d_int
            dy = y_d - y_d_int

            img_mat[y_d_int,x_d_int] = (1-dx)*(1-dy)*img[y,x] + dx*(1-dy)*img[y,x+1] + \
                               (1-dx)*dy*img[y+1,x] + dx*dy*img[y+1,x+1]
            img_mat[y_d_int,x_d_int+1] = (1-dx)*(1-dy)*img[y,x+1] + dx*(1-dy)*img[y,x+2] + \
                                    (1-dx)*dy*img[y+1,x+1] + dx*dy*img[y+1,x+2]
            img_mat[y_d_int+1,x_d_int] = (1-dx)*(1-dy)*img[y+1,x] + dx*(1-dy)*img[y+1,x+1] + \
                                    (1-dx)*dy*img[y+2,x] + dx*dy*img[y+2,x+1]
            img_mat[y_d_int+1,x_d_int+1] = (1-dx)*(1-dy)*img[y+1,x+1] + dx*(1-dy)*img[y+1,x+2] + \
                                    (1-dx)*dy*img[y+2,x+1] + dx*dy*img[y+2,x+2]


    return img_mat

Bilinear Interpolation output image

Method 2: Interpolation with scipy.ndimage.map_coordinates

In this method, I used scipy.ndimage.map_coordinates to apply radial distortion. I expected to get pincushion distortion, but the output showed barrel distortion. Here's the code :

def apply_distortion_scipy(img,W, k1, k2, k3):
    H = W
    x, y = np.meshgrid(np.float32(np.arange(W)), np.float32(np.arange(H)))  # meshgrid for interpolation mapping
    dis_center_x = W / 2
    dis_center_y = W / 2

    x_norm = (x - dis_center_x) / (W)
    y_norm = (y - dis_center_y) / (H)
    r = x_norm * x_norm + y_norm * y_norm

    coeff = (k1 + k2 * r + k3 * (r ** 2))
    x_d = coeff * (x - W / 2) + W / 2
    y_d = coeff * (y - H / 2) + H / 2

    # img = np.array(img)
    img_mat = scipy.ndimage.map_coordinates(img, [y_d.ravel(), x_d.ravel()])
    img_mat.resize(img.shape)


    return img_mat

scipy.ndimage.map_coordinates output image

I have looked through the API Reference of scipy.ndimage.map_coordinates, but there are no specific implementation codes. Can anyone tell me what could be wrong with my issue?

Solution

The issue is that coordinates argument of scipy.ndimage.map_coordinates applies the inverse coordinate map.
The result is the "inverted distortion".

The documentation is not very clear:

The array of coordinates is used to find, for each point in the output, the corresponding coordinates in the input.

The result is equivalent to the following implementation:

def apply_undistort_jit(img, W, img_mat, k1, k2, k3):
    H = W
    dis_center_x = W/2
    dis_center_y = H/2

    for x in range(W-2):
        for y in range(H-2):
            x_norm = (x - dis_center_x)/ (W)
            y_norm = (y - dis_center_y)/ (H)
            r = x_norm * x_norm + y_norm * y_norm

            coeff = (k1 + k2 * r + k3 * (r ** 2))
            x_d = coeff * (x - W/2) + W/2
            y_d = coeff * (y - H/2) + H/2
            if x_d >= W-1 or y_d >= H-1 or x_d<0 or y_d <0:
                continue

            x_d_int = math.floor(x_d)
            y_d_int = math.floor(y_d)

            dx = x_d - x_d_int
            dy = y_d - y_d_int

            #img_mat[y_d_int,x_d_int] = (1-dx)*(1-dy)*img[y,x] + dx*(1-dy)*img[y,x+1] + (1-dx)*dy*img[y+1,x] + dx*dy*img[y+1,x+1]
            #img_mat[y_d_int,x_d_int+1] = (1-dx)*(1-dy)*img[y,x+1] + dx*(1-dy)*img[y,x+2] + (1-dx)*dy*img[y+1,x+1] + dx*dy*img[y+1,x+2]
            #img_mat[y_d_int+1,x_d_int] = (1-dx)*(1-dy)*img[y+1,x] + dx*(1-dy)*img[y+1,x+1] + (1-dx)*dy*img[y+2,x] + dx*dy*img[y+2,x+1]
            #img_mat[y_d_int+1,x_d_int+1] = (1-dx)*(1-dy)*img[y+1,x+1] + dx*(1-dy)*img[y+1,x+2] + (1-dx)*dy*img[y+2,x+1] + dx*dy*img[y+2,x+2]
            img_mat[y, x] = (1-dx)*(1-dy)*img[y_d_int,x_d_int] + dx*(1-dy)*img[y_d_int,x_d_int+1] + (1-dx)*dy*img[y_d_int+1,x_d_int] + dx*dy*img[y_d_int+1,x_d_int+1]

Note that, all the geometrical transformation functions applies the same principe:

Iterate over the destination coordinates.
For each destination coordinate, look for the source coordinate.

The input to the geometrical transformation is the "backward map" (destination to source) and not the "forward map" (source to destination).

Note:
The apply_distortion_jit implementation incorrect.
Filling 4 pixels as implemented, causes visible artifacts:

img_mat[y_d_int,x_d_int] = (1-dx)*(1-dy)*img[y,x] + dx*(1-dy)*img[y,x+1] + (1-dx)*dy*img[y+1,x] + dx*dy*img[y+1,x+1]
img_mat[y_d_int,x_d_int+1] = (1-dx)*(1-dy)*img[y,x+1] + dx*(1-dy)*img[y,x+2] + (1-dx)*dy*img[y+1,x+1] + dx*dy*img[y+1,x+2]
img_mat[y_d_int+1,x_d_int] = (1-dx)*(1-dy)*img[y+1,x] + dx*(1-dy)*img[y+1,x+1] + (1-dx)*dy*img[y+2,x] + dx*dy*img[y+2,x+1]
img_mat[y_d_int+1,x_d_int+1] = (1-dx)*(1-dy)*img[y+1,x+1] + dx*(1-dy)*img[y+1,x+2] + (1-dx)*dy*img[y+2,x+1] + dx*dy*img[y+2,x+2]

The correct implementation (when having the inverse distortion) is:

img_mat[y, x] = (1-dx)*(1-dy)*img[y_d_int,x_d_int] + dx*(1-dy)*img[y_d_int,x_d_int+1] + (1-dx)*dy*img[y_d_int+1,x_d_int] + dx*dy*img[y_d_int+1,x_d_int+1]

I am going to use OpenCV instead of SciPy, because there are existing solutions using OpenCV.

According to the following answer, "OpenCV doesn't provide distort function for image".

The prefered solution is inverting the maps:

Invert the "remap" maps.
Apply cv2.remap (or scipy.ndimage.map_coordinates) using the inverted maps.

The following post shows ways for inverting the maps.

Code sample for "forward distortion":

import numpy as np
import math
import cv2

# https://stackoverflow.com/a/68706787/4926757
def invert_map(F):
    sh = (F.shape[1], F.shape[0])
    I = np.zeros_like(F)
    I[:,:,1], I[:,:,0] = np.indices(sh)
    P = np.copy(I)
    for i in range(10):
        P += I - cv2.remap(F, P, None, interpolation=cv2.INTER_LINEAR)
    return P

k_1 = 1
k_2 = 0.4
k_3 = 0.2

img = cv2.imread('input.jpg', cv2.IMREAD_GRAYSCALE)

dist_coeffs = np.array([k_1, k_2, 0, 0, k_3])  # (k1, k2, p1, p2, k3)
camera_matrix = np.eye(3)
camera_matrix[0, 0] = img.shape[1]
camera_matrix[1, 1] = img.shape[0]
camera_matrix[0, 2] = (img.shape[1]-1)/2
camera_matrix[1, 2] = (img.shape[0]-1)/2
new_camera_matrix = camera_matrix.copy()

# Compute "Undistort" maps:  
mapxy, _ = cv2.initUndistortRectifyMap(camera_matrix, dist_coeffs, None, new_camera_matrix, (img.shape[1], img.shape[0]), m1type=cv2.CV_32FC2)

# Invert the maps    
inv_mapxy = invert_map(mapxy)
inv_mapx = inv_mapxy[:, :, 0]
inv_mapy = inv_mapxy[:, :, 1]

# Use the inverted maps
out_img = cv2.remap(img, inv_mapx, inv_mapy, cv2.INTER_LINEAR)

cv2.imwrite('out_img.jpg', out_img)

Sample output (the input is your "flattened image" with the artifacts):