基于双线性插值的OpenCL仿射变换核Python重实现问题排查

我有一个实现带双线性插值的仿射变换的OpenCL核，想把它重写成Python代码。我已经完成了Python实现，但运行结果和期望不符，找不到问题所在。

OpenCL核代码 #

#ifndef SAMPLER_FILTER
#define SAMPLER_FILTER CLK_FILTER_LINEAR
#endif

#ifndef SAMPLER_ADDRESS
#define SAMPLER_ADDRESS CLK_ADDRESS_CLAMP
#endif

__kernel void affine_transform_2d_interpolate(
    IMAGE_input_TYPE input,
    IMAGE_output_TYPE output,
    IMAGE_mat_TYPE mat)
{

  const sampler_t sampler = CLK_NORMALIZED_COORDS_TRUE|
      SAMPLER_ADDRESS | SAMPLER_FILTER;

  uint i = get_global_id(0);
  uint j = get_global_id(1);

  uint Nx = GET_IMAGE_WIDTH(input);
  uint Ny = GET_IMAGE_HEIGHT(input);

  float x = i+0.5f;
  float y = j+0.5f;

  float y2 = (mat[3]*x+mat[4]*y+mat[5]);
  float x2 = (mat[0]*x+mat[1]*y+mat[2]);

  float2 coord_norm = (float2)(x2/Nx,y2/Ny);

  float pix = (float)(READ_input_IMAGE(input, sampler, coord_norm).x);
  int2 pos = (int2){i, j};

  WRITE_output_IMAGE(output, pos, CONVERT_output_PIXEL_TYPE(pix));
  
}

期望输出 #

我的Python重实现代码 #

import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imshow

def bilinear_interpolation(input_matrix, output_size):
    input_height, input_width = input_matrix.shape
    output_height, output_width = output_size

    x_ratio = (input_width - 1) / float(output_width - 1)
    y_ratio = (input_height - 1) / float(output_height - 1)

    output_matrix = np.zeros(output_size)

    input_height = input_height - 1
    input_width = input_width - 1 
    for y in range(output_height):
        for x in range(output_width):
            mat = [0.2, 0, 0, 0, 0.2, 0]

            x_input = mat[0]*x+mat[1]*y+mat[2]
            y_input = mat[3]*x+mat[4]*y+mat[5]
            
            # x_floor = int(np.floor(x_input))
            # y_floor = int(np.floor(y_input))
            x_floor = int(x_input)
            y_floor = int(y_input)

            x_weight = x_input - x_floor
            y_weight = y_input - y_floor
            
            top_left = input_matrix[x_floor, y_floor]
            top_right = input_matrix[(x_floor-1) if (x_floor+1)>input_width else (x_floor+1), (y_floor)]
            bottom_left = input_matrix[(x_floor), (y_floor-1) if (y_floor+1)>input_height else (y_floor+1)]
            bottom_right = input_matrix[(x_floor-1) if (x_floor+1)>input_width else (x_floor+1),
                                        (y_floor-1) if (y_floor+1)>input_height else (y_floor+1)]

            interpolated_value = (
                top_left * (1 - x_weight) * (1 - y_weight) +
                top_right * x_weight * (1 - y_weight) +
                bottom_left * (1 - x_weight) * y_weight +
                bottom_right * x_weight * y_weight
            )

            output_matrix[y, x] = interpolated_value

    return output_matrix

# Input 2x2 matrix
input_matrix = np.array([[1, 2], [3, 4]])

# Output 10x10 matrix
output_size = (10, 10)

# Perform bilinear interpolation
output_matrix = bilinear_interpolation(input_matrix, output_size)

imshow(output_matrix)

当前输出 #

问题排查与修复 #

对比OpenCL代码，Python实现存在以下核心错误：

1. 坐标基准错误：OpenCL使用像素中心坐标（i+0.5f、j+0.5f），Python代码直接用像素左上角索引值x、y，导致仿射变换坐标起点偏差。
2. 坐标映射错误：OpenCL将变换后的坐标x2、y2除以输入图像宽高得到归一化坐标，再由采样器映射到实际像素；Python代码直接用变换后的小数坐标计算，未对应输入图像的实际尺寸范围。
3. 数组索引顺序错误：NumPy数组索引是[行, 列]（对应图像y、x），OpenCL中图像宽是x方向、高是y方向，Python代码索引顺序完全颠倒。
4. 边界处理逻辑错误：OpenCL的CLK_ADDRESS_CLAMP是超出边界取最近边界像素，Python代码的条件判断和索引计算逻辑混乱，导致邻域像素选取错误。
5. 插值权重计算错误：基于错误坐标计算的权重无法得到正确插值结果。

修正后的Python代码 #

import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imshow

def affine_bilinear_interpolate(input_matrix, output_size, mat):
    input_height, input_width = input_matrix.shape
    output_height, output_width = output_size
    output_matrix = np.zeros(output_size, dtype=np.float32)
    
    for j in range(output_height):
        for i in range(output_width):
            # 对齐OpenCL的像素中心坐标
            x = i + 0.5
            y = j + 0.5
            
            # 应用仿射变换
            x2 = mat[0] * x + mat[1] * y + mat[2]
            y2 = mat[3] * x + mat[4] * y + mat[5]
            
            # 映射到输入图像的实际像素范围，并实现CLK_ADDRESS_CLAMP边界处理
            input_x = np.clip(x2, 0, input_width - 1)
            input_y = np.clip(y2, 0, input_height - 1)
            
            # 计算四个邻域像素坐标
            x_floor = np.floor(input_x).astype(int)
            y_floor = np.floor(input_y).astype(int)
            x_ceil = min(x_floor + 1, input_width - 1)
            y_ceil = min(y_floor + 1, input_height - 1)
            
            # 按NumPy索引顺序[y, x]获取邻域像素值
            top_left = input_matrix[y_floor, x_floor]
            top_right = input_matrix[y_floor, x_ceil]
            bottom_left = input_matrix[y_ceil, x_floor]
            bottom_right = input_matrix[y_ceil, x_ceil]
            
            # 计算插值权重
            x_weight = input_x - x_floor
            y_weight = input_y - y_floor
            
            # 双线性插值计算
            interpolated_value = (
                top_left * (1 - x_weight) * (1 - y_weight) +
                top_right * x_weight * (1 - y_weight) +
                bottom_left * (1 - x_weight) * y_weight +
                bottom_right * x_weight * y_weight
            )
            
            output_matrix[j, i] = interpolated_value
    
    return output_matrix

# 输入2x2矩阵
input_matrix = np.array([[1, 2], [3, 4]], dtype=np.float32)
# 输出10x10矩阵
output_size = (10, 10)
# 仿射变换矩阵（与OpenCL中的mat一致）
mat = [0.2, 0, 0, 0, 0.2, 0]

# 执行变换
output_matrix = affine_bilinear_interpolate(input_matrix, output_size, mat)

imshow(output_matrix, cmap='gray')
plt.show()

修复说明 #

• 修正坐标基准：使用像素中心坐标，和OpenCL逻辑完全对齐。
• 对齐坐标映射：直接使用变换后的坐标，通过np.clip实现CLK_ADDRESS_CLAMP边界处理。
• 调整索引顺序：按NumPy的[y, x]索引访问像素，匹配图像行列对应关系。
• 简化邻域选取：直接计算上下左右邻域，超出边界时取边界值。
• 提取变换矩阵为参数：让函数更灵活，与OpenCL接口逻辑一致。

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