Scipy.ndimage是否存在与map_coordinates反向的推送式变形函数？

scipy.ndimage的map_coordinates是拉取式的变形操作：它接收源体数据和目标空间映射（指定目标空间每个体素对应的源空间位置），从源数据中拉取对应位置的值生成结果。

你需要的是反向的推送式变形：输入源体数据和变形向量映射（每个源值对应的移动位置向量），把源数据的每个元素按向量推送到目标位置，同时支持插值，用来生成中间阶段数据。

先回顾map_coordinates的用法示例 #

from scipy import ndimage
import numpy as np
a = np.arange(4.).reshape((2, 2))
print(a)

输出：

array([[ 0.,  1.],
       [ 2.,  3.]])

print(ndimage.map_coordinates(a, [[0.5, 1], [0.5, 1]], order=1))

输出：

array([ 1.5,  3])

你期望的推送式变形示例 #

from scipy import ndimage
import numpy as np
a = np.arange(4.).reshape((2, 2))
print(a)

输出：

array([[ 0.,  1.],
       [ 2.,  3.]])

deformation = np.array([[[0.5,0.5],[1,1]],
                         [[0,0.5],[1,0.5]]])
apply_deformation_with_interpolation(a, deformation)

解答 #

Scipy.ndimage本身没有直接提供这种“推送”式的变形函数，但可以通过两种思路实现等效效果：

思路1：反向映射+map_coordinates #

推送式变形的核心逻辑是：源位置s的元素会移动到t = s + 变形向量[s]。反过来，目标位置t的取值需要找到所有会移动到这里的源位置并加权求和。我们可以构建反向映射——对每个目标坐标t，计算对应的源坐标s（满足s + 变形向量[s] = t），然后用map_coordinates从源数据拉取s处的值。

如果变形幅度较小，可以近似用s ≈ t - 变形向量[t]来简化计算，直接基于目标坐标反向推导源坐标。

思路2：手动实现Splatting（喷洒）操作 #

另一种方法是直接模拟“推送”过程：遍历每个源位置，将对应的值按变形向量推送到目标位置，同时用插值（比如双线性插值）把值分配到目标位置周围的体素上。以下是一个2D场景的简单实现：

import numpy as np

def apply_deformation_with_interpolation(source, deformation):
    target = np.zeros_like(source)
    source_shape = source.shape
    # 展开所有源位置的坐标和对应值
    source_coords = np.indices(source_shape).reshape(len(source_shape), -1)
    target_coords = source_coords + deformation.reshape(len(source_shape), -1)
    values = source.reshape(-1)
    
    # 双线性插值分配值到目标体素
    if source.ndim == 2:
        for (sx, sy), val in zip(source_coords.T, values):
            tx, ty = target_coords[:, sx * source_shape[1] + sy]
            # 获取目标坐标周围的整数网格点
            x0, x1 = int(np.floor(tx)), int(np.ceil(tx))
            y0, y1 = int(np.floor(ty)), int(np.ceil(ty))
            # 计算各网格点的权重
            wx1, wx0 = tx - x0, x1 - tx
            wy1, wy0 = ty - y0, y1 - ty
            
            # 将值加权分配到有效范围内的网格点
            for x in [x0, x1]:
                for y in [y0, y1]:
                    if 0 <= x < source_shape[0] and 0 <= y < source_shape[1]:
                        weight = wx0 if x == x0 else wx1
                        weight *= wy0 if y == y0 else wy1
                        target[x, y] += val * weight
    return target

# 测试代码
a = np.arange(4.).reshape((2, 2))
print("源数据：")
print(a)
deformation = np.array([[[0.5,0.5],[1,1]],
                        [[0,0.5],[1,0.5]]])
result = apply_deformation_with_interpolation(a, deformation)
print("\n推送变形结果：")
print(result)

如果需要处理大变形、高维数据或者更高效的实现，可以考虑使用SimpleITK、Elastix这类专门的影像处理库，它们提供了更完善的推送式变形支持。

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