如何使用Numpy获取二维矩阵的直线采样系数（最近邻法）？

是否有现成的Numpy/Scipy函数可实现此功能（例如Bresenham算法）？我无法确定以下的简易版本能否适用于所有场景。

import numpy as np

a = np.arange(36).reshape((6, 6))
# [[ 0  1  2  3  4  5]
#  [ 6  7  8  9 10 11]
#  [12 13 14 15 16 17]
#  [18 19 20 21 22 23]
#  [24 25 26 27 28 29]
#  [30 31 32 33 34 35]]

def sample_over_line(a, i0, j0, i1, j1):
    di = i1 - i0
    dj = j1 - j0
    if abs(di) > abs(dj):
        if i0 > i1:
            i0, j0, i1, j1 = i1, j1, i0, j0
        b = [a[i, round((j1 - j0) / (i1 - i0) * (i - i0) + j0)] for i in range(i0, i1 + 1)]
    else:
        if j0 > j1:
            i0, j0, i1, j1 = i1, j1, i0, j0
        b = [a[round((i1 - i0) / (j1 - j0) * (j - j0) + i0), j] for j in range(j0, j1 + 1)]
    return b

print(sample_over_line(a=a, i0=1, j0=0, i1=3, j1=5))  # 从[1,0]到[3,5]的直线采样结果
# [6, 7, 14, 15, 22, 23]

print(sample_over_line(a=a, i0=1, j0=2, i1=5, j1=3))  # 从[1,2]到[5,3]的直线采样结果
# [8, 14, 20, 27, 33]

print(sample_over_line(a=a, i0=5, j0=3, i1=1, j1=2)) # 从[5,3]到[1,2]的直线采样结果
# [8, 14, 20, 27, 33]

解决方案 #

1. 用Scipy工具快速实现

Scipy的ndimage.map_coordinates支持最近邻插值，能直接满足需求，且无需手动处理边界情况：

import numpy as np
from scipy.ndimage import map_coordinates

a = np.arange(36).reshape((6, 6))

def sample_line_scipy(a, i0, j0, i1, j1):
    # 生成直线上的采样点坐标
    num_points = max(abs(i1-i0), abs(j1-j0)) + 1
    i_coords = np.linspace(i0, i1, num_points)
    j_coords = np.linspace(j0, j1, num_points)
    # 最近邻采样（order=0）
    return map_coordinates(a, [i_coords, j_coords], order=0).astype(int)

# 测试
print(sample_line_scipy(a, 1, 0, 3, 5))       # [ 6  7 14 15 22 23]
print(sample_line_scipy(a, 1, 2, 5, 3))       # [ 8 14 20 27 33]
print(sample_line_scipy(a, 5, 3, 1, 2))       # [33 27 20 14  8]

2. 手动实现Bresenham算法（无依赖）

如果不想依赖Scipy，标准Bresenham算法能高效生成直线上的整数坐标，覆盖所有场景（包括垂直、水平直线）：

import numpy as np

def bresenham(i0, j0, i1, j1):
    points = []
    di = abs(i1 - i0)
    dj = abs(j1 - j0)
    i, j = i0, j0
    si = 1 if i1 > i0 else -1
    sj = 1 if j1 > j0 else -1
    
    if di > dj:
        err = di / 2.0
        while i != i1 + si:
            points.append((i, j))
            err -= dj
            if err < 0:
                j += sj
                err += di
            i += si
    else:
        err = dj / 2.0
        while j != j1 + sj:
            points.append((i, j))
            err -= di
            if err < 0:
                i += si
                err += dj
            j += sj
    return points

def sample_line_bresenham(a, i0, j0, i1, j1):
    coords = bresenham(i0, j0, i1, j1)
    return [a[i, j] for i, j in coords]

# 测试
a = np.arange(36).reshape((6, 6))
print(sample_line_bresenham(a, 1, 0, 3, 5))   # [6, 7, 14, 15, 22, 23]
print(sample_line_bresenham(a, 1, 2, 5, 3))   # [8, 14, 20, 27, 33]
print(sample_line_bresenham(a, 5, 3, 1, 2))   # [33, 27, 20, 14, 8]

3. 现有代码的问题

你当前的实现存在几个隐患：

• 垂直/水平直线会触发除零错误；
• 未覆盖斜率为0或无穷大的边界场景；
• 列表推导式效率低于Numpy向量操作，大矩阵场景下表现差。

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