Как сделать 2D-разрез через 3D-изображение?

Question

У меня есть 3D-массив с некоторыми данными (растровое 3D-изображение). Я хотел бы получить 2D-разрез через этот массив, используя некоторую подходящую интерполяцию (предпочтительно линейную - в данном случае, вероятно, "трилинейную"). Плоскость среза может быть описана, однако удобно, например, с использованием вектора нормали и расстояния.

Если разрез параллелен одной из осей, это тривиально, просто нарежьте трехмерный массив (с кусочком индексного среза). Но если разрез не параллелен оси, я не вижу хорошего способа начать с этой проблемой. Единственное, что приходит на ум, - это вращать трехмерный массив (возможно, используя композицию двумерных вращений) так, чтобы разрез проходил параллельно оси, но это кажется ужасно неэффективным.

Я работаю в Python с numpy, ndimage и skimage. Можно предположить, что доступны любые другие модули Python.

Answer 1

На самом деле это не проверялось, но это действительно производит какое-то изображение. По предложению @Daniel Forsman.

import numpy as np
from scipy.interpolate import RegularGridInterpolator

# stack coordinates
z0,z1,z2 = 20, 20, 20
zz0,zz1,zz2 = np.linspace(0, 1, z0), np.linspace(0, 1, z1), np.linspace(0, 1, z2)

# fake stack data
d0,d1,d2 = np.ix_(0.5-np.abs(zz0-0.5), 0.5-np.abs(zz1-0.5), 0.5-np.abs(zz2-0.5))
data = np.minimum(np.minimum(d0, d1), d2)

# define picture (same coords as stack)
tl = np.array((0.1, -0.02, 0.3)) # top left corner
yo = np.array((-0.01, 0.1, 0.01))
yo /= np.sqrt((yo*yo).sum()) # y-axis unit
xo = np.array((0.1, 0, 0.1))
xo -= (xo*yo).sum() * yo # should be perpendicular now
xo /= np.sqrt((xo*xo).sum()) # x-axis unit

# build picture grid
nx,ny = 20j, 20j
ya, xa = np.ogrid[:1:ny, :1:nx]
grid = tl + ya[..., None] * yo + xa[..., None] * xo

picture = RegularGridInterpolator((zz0,zz1,zz2), data, bounds_error=False)(grid)