用 numpy 和/或 scipy 插入 3D 体积

Posted 2023-02-16

【中文标题】用 numpy 和/或 scipy 插入 3D 体积

【英文标题】：interpolate 3D volume with numpy and or scipy

【发布时间】：2014-03-17 03:52:30

【问题描述】：

我非常沮丧，因为几个小时后，我似乎无法在 python 中进行看似简单的 3D 插值。在 Matlab 中，我所要做的就是

Vi = interp3(x,y,z,V,xi,yi,zi)

使用 scipy 的 ndimage.map_coordinate 或其他 numpy 方法与此完全等价的是什么？

谢谢

【参考方案1】：

在 scipy 0.14 或更高版本中，有一个新函数 scipy.interpolate.RegularGridInterpolator 与 interp3 非常相似。

MATLAB 命令Vi = interp3(x,y,z,V,xi,yi,zi) 将转换为：

from numpy import array
from scipy.interpolate import RegularGridInterpolator as rgi
my_interpolating_function = rgi((x,y,z), V)
Vi = my_interpolating_function(array([xi,yi,zi]).T)

这是一个完整的例子，展示了两者；它将帮助您了解确切的差异...

MATLAB 代码：

x = linspace(1,4,11);
y = linspace(4,7,22);
z = linspace(7,9,33);
V = zeros(22,11,33);
for i=1:11
    for j=1:22
        for k=1:33
            V(j,i,k) = 100*x(i) + 10*y(j) + z(k);
        end
    end
end
xq = [2,3];
yq = [6,5];
zq = [8,7];
Vi = interp3(x,y,z,V,xq,yq,zq);

结果是Vi=[268 357]，这确实是(2,6,8) 和(3,5,7) 这两个点的值。

SCIPY 代码：

from scipy.interpolate import RegularGridInterpolator
from numpy import linspace, zeros, array
x = linspace(1,4,11)
y = linspace(4,7,22)
z = linspace(7,9,33)
V = zeros((11,22,33))
for i in range(11):
    for j in range(22):
        for k in range(33):
            V[i,j,k] = 100*x[i] + 10*y[j] + z[k]
fn = RegularGridInterpolator((x,y,z), V)
pts = array([[2,6,8],[3,5,7]])
print(fn(pts))

又是[268,357]。所以你会看到一些细微的差别：Scipy 使用 x,y,z 索引顺序，而 MATLAB 使用 y,x,z（奇怪）；在 Scipy 中，您在单独的步骤中定义一个函数，当您调用它时，坐标被分组为 (x1,y1,z1),(x2,y2,z2),... 而 matlab 使用 (x1,x2,.. .),(y1,y2,...),(z1,z2,...)。

除此之外，两者相似且同样易于使用。

【参考方案2】：

精确相当于 MATLAB 的 interp3 将使用 scipy 的 interpn 进行一次性插值：

import numpy as np
from scipy.interpolate import interpn

Vi = interpn((x,y,z), V, np.array([xi,yi,zi]).T)

MATLAB 和 scipy 的默认方法是线性插值，这可以通过 method 参数进行更改。请注意，interpn 仅支持 3 维及以上的线性和最近邻插值，这与支持三次和样条插值的 MATLAB 不同。

在同一网格上进行多次插值调用时，最好使用插值对象RegularGridInterpolator，如接受的答案above。 interpn 在内部使用 RegularGridInterpolator。

【参考方案3】：

基本上，ndimage.map_coordinates 在“索引”坐标（也称为“体素”或“像素”坐标）中工作。起初它的界面似乎有点笨拙，但它确实为您提供了很多的灵活性。

如果要指定类似于 matlab 的 interp3 的插值坐标，则需要将输入坐标转换为“索引”坐标。

还有一个额外的问题是map_coordinates 总是在输出中保留输入数组的 dtype。如果你插入一个整数数组，你会得到整数输出，这可能是也可能不是你想要的。对于下面的代码 sn-p，我假设您总是想要浮点输出。 （如果你不这样做，它实际上更简单。）

今晚晚些时候我会尝试添加更多解释（这是相当密集的代码）。

总而言之，我拥有的interp3 函数比您的确切目的可能需要的要复杂得多。但是，它应该或多或少地复制了我记得的 interp3 的行为（忽略 interp3(data, zoom_factor) 的“缩放”功能，scipy.ndimage.zoom 处理。）

import numpy as np
from scipy.ndimage import map_coordinates

def main():
    data = np.arange(5*4*3).reshape(5,4,3)

    x = np.linspace(5, 10, data.shape[0])
    y = np.linspace(10, 20, data.shape[1])
    z = np.linspace(-100, 0, data.shape[2])

    # Interpolate at a single point
    print interp3(x, y, z, data, 7.5, 13.2, -27)

    # Interpolate a region of the x-y plane at z=-25
    xi, yi = np.mgrid[6:8:10j, 13:18:10j]
    print interp3(x, y, z, data, xi, yi, -25 * np.ones_like(xi))

def interp3(x, y, z, v, xi, yi, zi, **kwargs):
    """Sample a 3D array "v" with pixel corner locations at "x","y","z" at the
    points in "xi", "yi", "zi" using linear interpolation. Additional kwargs
    are passed on to ``scipy.ndimage.map_coordinates``."""
    def index_coords(corner_locs, interp_locs):
        index = np.arange(len(corner_locs))
        if np.all(np.diff(corner_locs) < 0):
            corner_locs, index = corner_locs[::-1], index[::-1]
        return np.interp(interp_locs, corner_locs, index)

    orig_shape = np.asarray(xi).shape
    xi, yi, zi = np.atleast_1d(xi, yi, zi)
    for arr in [xi, yi, zi]:
        arr.shape = -1

    output = np.empty(xi.shape, dtype=float)
    coords = [index_coords(*item) for item in zip([x, y, z], [xi, yi, zi])]

    map_coordinates(v, coords, order=1, output=output, **kwargs)

    return output.reshape(orig_shape)

main()

【参考方案4】：

这个问题很老，但我认为它需要澄清一下，因为没有人指出请求的操作 (trilinear interpolation) 可以很容易地从头开始实现，同时持续节省计算时间（大约快 10 倍）w.r.t. scipy.interpolate的RegularGridInterpolator。

代码

import numpy as np
from itertools import product

def trilinear_interpolation(x_volume, y_volume, z_volume, volume, x_needed, y_needed, z_needed):
    """
    Trilinear interpolation (from Wikipedia)

    :param x_volume: x points of the volume grid 
    :type crack_type: list or numpy.ndarray
    :param y_volume: y points of the volume grid 
    :type crack_type: list or numpy.ndarray
    :param x_volume: z points of the volume grid 
    :type crack_type: list or numpy.ndarray
    :param volume:   volume
    :type crack_type: list or numpy.ndarray
    :param x_needed: desired x coordinate of volume
    :type crack_type: float
    :param y_needed: desired y coordinate of volume
    :type crack_type: float
    :param z_needed: desired z coordinate of volume
    :type crack_type: float

    :return volume_needed: desired value of the volume, i.e. volume(x_needed, y_needed, z_needed)
    :type volume_needed: float
    """
    # dimensinoal check
    if np.shape(volume) != (len(x_volume), len(y_volume), len(z_volume)):
        raise ValueError(f'dimension mismatch, volume must be a (len(x_volume), len(y_volume), len(z_volume)) list or numpy.ndarray')
    # check of the indices needed for the correct control volume definition
    i = searchsorted(x_volume, x_needed)
    j = searchsorted(y_volume, y_needed)
    k = searchsorted(z_volume, z_needed)
    # control volume definition
    control_volume_coordinates = np.array(
        [[x_volume[i - 1], y_volume[j - 1], z_volume[k - 1]], [x_volume[i], y_volume[j], z_volume[k]]])
    xd = (np.array([x_needed, y_needed, z_needed]) - control_volume_coordinates[0]) / (control_volume_coordinates[1] - control_volume_coordinates[0])
    # interpolation along x
    c2 = [[0, 0], [0, 0]]
    for m, n in product([0, 1], [0, 1]):
        c2[m][n] = volume[i - 1][j - 1 + m][k - 1 + n] * (1 - xd[0]) + volume[i][j - 1 + m][k - 1 + n] * xd[0]
    # interpolation along y
    c1 = [0, 0]
    c1[0] = c2[0][0] * (1 - xd[1]) + c2[1][0] * xd[1]
    c1[1] = c2[0][1] * (1 - xd[1]) + c2[1][1] * xd[1]
    # interpolation along z
    volume_needed = c1[0] * (1 - xd[2]) + c1[1] * xd[2]
    return volume_needed

def searchsorted(l, x):
    for i in l:
        if i >= x: break
    return l.index(i)


from scipy.interpolate import RegularGridInterpolator
def trilin_interp_regular_grid(x_volume, y_volume, z_volume, volume, x_needed, y_needed, z_needed):
    # dimensinoal check
    if np.shape(volume) != (len(x_volume), len(y_volume), len(z_volume)):
        raise ValueError(f'dimension mismatch, volume must be a (len(x_volume), len(y_volume), len(z_volume)) list or numpy.ndarray')
    # trilinear interpolation on a regular grid
    fn = RegularGridInterpolator((x_volume,y_volume,z_volume), volume)
    volume_needed = fn(np.array([x_needed, y_needed, z_needed]))
    return volume_needed

示例

import numpy as np
import time

the_volume = np.array(
[[[0.902, 0.985, 1.12, 1.267, 1.366],
[0.822, 0.871, 0.959, 1.064, 1.141],
[0.744, 0.77, 0.824, 0.897, 0.954],
[0.669, 0.682, 0.715, 0.765, 0.806],
[0.597, 0.607, 0.631, 0.667, 0.695]],
[[1.059, 1.09, 1.384, 1.682, 1.881],
[0.948, 0.951, 1.079, 1.188, 1.251],
[0.792, 0.832, 0.888, 0.940, 0.971],
[0.726, 0.733, 0.754, 0.777, 0.792],
[0.642, 0.656, 0.675, 0.691, 0.700]]])

x_needed = np.linspace(100, 1000, 2)
y_needed = np.linspace(0.3, 1, 3)
z_needed = np.linspace(0, 1, 3)
start = time.time()
manual_trilint = []
for x in x_needed:
    for y in y_needed:
        for z in z_needed:
            manual_trilint.append(trilinear_interpolation([100, 1000], [0.2, 0.4, 0.6, 0.8, 1], [0, 0.2, 0.5, 0.8, 1], the_volume, x, y, z))
end = time.time()
print(end - start)

start = time.time()
auto_trilint = []
for x in x_needed:
    for y in y_needed:
        for z in z_needed:
            auto_trilint.append(trilin_interp_regular_grid([100, 1000], [0.2, 0.4, 0.6, 0.8, 1], [0, 0.2, 0.5, 0.8, 1], the_volume, x, y, z))
end = time.time()
print(end - start)

【讨论】：

在插入大量点时，这实际上不是很省时。这可以从广播中受益，从而获得很大的加速