python - 图像处理中的径向畸变问题 - 使用桶形畸变而不是枕形畸变

Question

我在图像处理中遇到径向畸变问题。我正在尝试使用畸变系数 k_1 = 1 、 k_2 = 0.4 和 k_3 = 0.2 将径向畸变应用于图像。根据失真模型，我应该得到枕形失真作为输出。然而，当我尝试两种不同的方法时，我得到了两种不同的结果。

方法一:数组计算双线性插值

在这个方法中，我使用纯数组计算来应用径向畸变，输出是枕形畸变。这是我使用的代码:

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

方法2:使用 scipy.ndimage.map_coordinates 进行插值

在此方法中，我使用 scipy.ndimage.map_coordinates 来应用径向畸变。我预计会出现枕形失真，但输出显示桶形失真。代码如下:

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

我查阅了scipy.ndimage.map_coordinates的API Reference，但没有具体的实现代码。谁能告诉我我的问题可能出了什么问题？

Answer 1

问题在于 scipy.ndimage.map_coordinates 的 coefficients 参数应用了逆坐标图。
结果就是“反转失真”。

documentation不是很清楚:

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

结果相当于以下实现:

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]

请注意，所有几何变换函数都应用相同的原理:

• 迭代目标坐标。
• 对于每个目标坐标，查找源坐标。

几何变换的输入是“向后映射”(目的地到源)，而不是“向前映射”(源到目的地)。

---

注意:
apply_ Distortion_jit 实现不正确。
按照实现方式填充 4 个像素会导致可见的伪像:

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]  

---

我将使用 OpenCV 而不是 SciPy，因为已有使用 OpenCV 的解决方案。

根据following answer ，“OpenCV不提供图像扭曲功能”。

首选的解决方案是反转 map :

• 反转“重映射”映射。
• 使用倒置贴图应用cv2.remap(或scipy.ndimage.map_coordinates)。

following post显示反转 map 的方法。

---

“前向失真”的代码示例:

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)

---

示例输出(输入是带有伪影的“扁平图像”):