图像处理径向畸变异常：Scipy输出桶形而非预期枕形畸变

我在图像处理中尝试应用径向畸变，使用畸变系数k₁=1、k₂=0.4、k₃=0.2，根据畸变模型应得到枕形畸变，但两种实现方法结果不同：

方法1：纯数组计算双线性插值 #

通过纯数组计算实现径向畸变，输出为预期的枕形畸变，代码如下：

import math

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

输出为预期的枕形畸变图像。

方法2：使用scipy.ndimage.map_coordinates插值 #

使用该方法实现径向畸变，预期得到枕形畸变，但输出为桶形畸变，代码如下：

import numpy as np
import scipy.ndimage

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_mat = scipy.ndimage.map_coordinates(img, [y_d.ravel(), x_d.ravel()])
    img_mat.resize(img.shape)

    return img_mat

输出为桶形畸变图像。

问题原因 #

核心差异在于映射方向的逻辑相反：

• 方法1采用正向映射：遍历原始图像的每个像素，计算它在畸变后图像中的目标位置，将原始像素值插值填充到该位置。当coeff>1时，原始边缘像素会被映射到更外侧的位置，最终形成枕形畸变。
• scipy.ndimage.map_coordinates采用反向映射逻辑：它需要接收的是「畸变后图像的每个像素，对应到原始图像中的坐标位置」，而你当前代码传入的是「原始图像像素映射到畸变后图像的坐标」，相当于把映射方向完全搞反了，最终得到了与预期相反的桶形畸变。

修正思路 #

需要计算反向映射关系：对于畸变后图像的每个位置(x_d, y_d)，找到它在原始图像中对应的坐标(x, y)，具体步骤如下：

1. 生成畸变后图像的网格坐标(x_d, y_d)（即当前代码中的x, y）
2. 计算该位置的归一化坐标：x_norm_d = (x_d - dis_center_x)/W，y_norm_d = (y_d - dis_center_y)/H，进而得到r_d = x_norm_d² + y_norm_d²
3. 根据径向畸变模型，求解原始归一化半径r：满足 r_d = r*(k1 + k2*r² + k3*r⁴)（可通过数值方法求解，比如牛顿迭代）
4. 计算原始坐标：

   if r_d != 0:
       x = dis_center_x + r * (x_d - dis_center_x) / np.sqrt(r_d)
       y = dis_center_y + r * (y_d - dis_center_y) / np.sqrt(r_d)
   else:
       x, y = dis_center_x, dis_center_y
   

5. 将计算得到的原始坐标[y.ravel(), x.ravel()]传入map_coordinates，即可得到正确的枕形畸变结果。

内容的提问来源于stack exchange，提问作者Elanie