Python：快速离散插值

Question

给定一个表示离散坐标变换函数的 4D 数组

arr[x_in, y_in, z_in] = [x_out, y_out, z_out]

我想将 arr 插入到具有更多元素的网格中（假设 arr 中的样本最初是从高元素立方体的规则间隔网格中抽取的）。我已经尝试了 scipy 中的 RegularGridInterpolator，但是这相当慢：

import numpy as np
from scipy.interpolate import RegularGridInterpolator
from time import time

target_size   = 32
reduced_size  = 5

small_shape = (reduced_size,reduced_size,reduced_size,3)
cube_small  = np.random.randint(target_size, size=small_shape, dtype=np.uint8)

igrid = 3*[np.linspace(0, target_size-1, reduced_size)]
large_shape = (target_size, target_size, target_size,3)
cube_large  = np.empty(large_shape)

t0 = time()
interpol = RegularGridInterpolator(igrid, cube_small)
for x in np.arange(target_size):
    for y in np.arange(target_size):
        for z in np.arange(target_size):
            cube_large[x,y,z] = interpol([x,y,z])
print(time()-t0)

有没有想到更适合该任务的算法？也许有一些东西可以利用这样一个事实，在这种情况下，我只对每个点的离散整数值感兴趣。

Answer 1

对于网格的生成我真的不能说，因为我不完全确定你正在尝试做什么。但是就提高效率而言，使用三重 for 循环而不是使用广播会大大降低您的代码速度

import itertools
import numpy as np
from scipy.interpolate import RegularGridInterpolator
from time import time

target_size   = 32
reduced_size  = 5

small_shape = (reduced_size,reduced_size,reduced_size,3)
cube_small  = np.random.randint(target_size, size=small_shape, dtype=np.uint8)


igrid = 3*[np.linspace(0, target_size-1, reduced_size)]


large_shape = (target_size, target_size, target_size,3)
cube_large  = np.empty(large_shape)

def original_method():
    t0 = time()

    interpol = RegularGridInterpolator(igrid, cube_small)
    for x in np.arange(target_size):
        for y in np.arange(target_size):
            for z in np.arange(target_size):
                cube_large[x,y,z] = interpol([x,y,z])

    print('ORIGINAL METHOD: ', time()-t0)
    return cube_large

def improved_method():
    t1 = time()
    interpol = RegularGridInterpolator(igrid, cube_small)

    arr = np.arange(target_size)
    grid = np.array(list(itertools.product(arr, repeat=3)))
    cube_large = interpol(grid).reshape(target_size, target_size, target_size, 3)

    print('IMPROVED METHOD:', time() - t1)

    return cube_large


c1 = original_method()
c2 = improved_method()

print('Is the result the same?  ', np.all(c1 == c2))

输出

ORIGINAL METHOD:  6.9040000438690186
IMPROVED METHOD: 0.026999950408935547
Is the result the same?   True

Answer 2

确实按照建议通过了整个网格似乎更有效。我正在添加一个 numpy 版本，它似乎比 itertools:

快 2 倍

target_size = 256

... # Same as above

def do_new():
    t0 = time()
    grid = np.meshgrid(*3*[np.arange(target_size)], indexing='ij')
    grid = grid[::-1]
    grid = np.stack(grid).T
    ret  = interpol(grid)
    print(f"{time()-t0}")
    return ret

c1 = improved_method()
c2 = do_new()
print((c1==c2).all())

输出：

IMPROVED METHOD: 13.865211009979248
6.268443584442139
True

Python：快速离散插值

Python: Fast discrete interpolation

python

interpolation

numpy

scipy