1

I have a 3D array with some data (a raster 3D image). I would like to get a 2D cut through that array, using some suitable interpolation (preferably linear - that's probably "trilinear" in this case). The plane of the cut can be described however is convenient, for example using a normal vector and distance.

If the cut is parallel to one of the axes, this is trivial, just slice the 3D array (with numpy index slice). But if the cut is not parallel to an axis, I don't see a good way to get started with that problem. The only thing that comes to mind is to rotate the 3D array (probably using a composition of 2D rotations) so that the cut is parallel to an axis, but that seems terribly inefficient.

I am working in python with numpy, ndimage and skimage. Any other python modules may be assumed to be available.

Alex I
  • 19,689
  • 9
  • 86
  • 158
  • 1
    You could rotate your planar coordinates and use [`scipy.interpolate.RegularGridInterpolator`](https://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.interpolate.RegularGridInterpolator.html) – Daniel F Jan 30 '17 at 07:58
  • @DanielForsman Thanks a lot, that sounds like it would do the job. You can post some example code as an answer if you'd like, I'll accept it – Alex I Jan 30 '17 at 08:49
  • Don't really know enough to make a complete answer, it was more just google-fu. – Daniel F Jan 30 '17 at 09:03

1 Answers1

2

Didn't really test this but it does produce an image of sorts. Based on @Daniel Forsman's suggestion.

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)
Paul Panzer
  • 51,835
  • 3
  • 54
  • 99