Repository navigation
[Doc]: Gallery example showing 3D slice planes #28239
Description
@timhoffm suggested in #3919 (comment) that we make a gallery example for plotting X Y Z planar slices through a volume.
His result and its code:
import matplotlib.pyplot as plt
import numpy as np
def plot_quadrants(ax, array, fixed_coord, cmap):
"""For a given 3d *array* plot a plane with *fixed_coord*, using four individual quadrants."""
nx, ny, nz = array.shape
index = {
'x': (nx//2, slice(None), slice(None)),
'y': (slice(None), ny//2, slice(None)),
'z': (slice(None), slice(None), nz//2),
}[fixed_coord]
plane_data = array[index]
n0, n1 = plane_data.shape
quadrants = [
plane_data[:n0//2, :n1//2],
plane_data[:n0//2, n1//2:],
plane_data[n0//2:, :n1//2],
plane_data[n0//2:, n1//2:]
]
min_val = array.min()
max_val = array.max()
cmap = plt.get_cmap(cmap)
for i, quadrant in enumerate(quadrants):
facecolors = cmap((quadrant - min_val) / (max_val-min_val))
if fixed_coord == 'x':
Y, Z = np.mgrid[0:ny//2, 0:nz//2]
X = nx//2 * np.ones_like(Y)
Y_offset = (i // 2) * ny//2
Z_offset = (i % 2) * nz//2
ax.plot_surface(X, Y + Y_offset, Z + Z_offset, rstride=1, cstride=1, facecolors=facecolors, shade=False)
elif fixed_coord == 'y':
X, Z = np.mgrid[0:nx//2, 0:nz//2]
Y = ny//2 * np.ones_like(X)
X_offset = (i // 2) * nx//2
Z_offset = (i % 2) * nz//2
ax.plot_surface(X + X_offset, Y, Z + Z_offset, rstride=1, cstride=1, facecolors=facecolors, shade=False)
elif fixed_coord == 'z':
X, Y = np.mgrid[0:nx//2, 0:ny//2]
Z = nz//2 * np.ones_like(X)
X_offset = (i // 2) * nx//2
Y_offset = (i % 2) * ny//2
ax.plot_surface(X + X_offset, Y + Y_offset, Z, rstride=1, cstride=1, facecolors=facecolors, shade=False)
def figure_3D_array_slices(array, cmap=None):
"""Plot a 3d array using three intersecting centered planes."""
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.set_box_aspect(array.shape)
plot_quadrants(ax, array, 'x', cmap=cmap)
plot_quadrants(ax, array, 'y', cmap=cmap)
plot_quadrants(ax, array, 'z', cmap=cmap)
return fig, ax
nx, ny, nz = 70, 100, 50
r_square = (np.mgrid[-1:1:1j*nx, -1:1:1j*ny, -1:1:1j*nz]**2).sum(0)
figure_3D_array_slices(r_square, cmap='viridis_r')
plt.show()Reacted by Tim Hoffmann
Activity
story645 commentedon May 17, 2024
Um can it be split up into three subplots, one for each axes?
Just thinking through that while pretty, it might be a bit hard to follow from a "what do I do to make this slice for one axes"?
timhoffm commentedon May 17, 2024
No. The point here is to illustrate multiple intersecting planes in one subplot, which requires the quadrant subdivision for correct drawing order.
This is somewhat inherently complex (with current API (?)) due to the involved indexing. Note also that indexing will get even worse if you don’t split in the center (which I thus have refained from to generalize to). I rather see this as an illustration/proof of concept how you can work around limitations of the 3d rendering.
Separate planes (well, one would be enough then) is a simpler example that should be done separately. Also note that drawing a single color-coded plane is currently already quite a hassle because you have to use
plot_surface. Somequadmeshlike function to draw axes-aligned planes from 2D arrays would be reasonable. But that’s even a bigger topic.story645 commentedon May 17, 2024
My concern more is that in this case the complexity/way that it's written may make it hard to unravel and tweak for someone who wanted to make something similar but not identical.
I wonder if this would work better as a tutorial where these stages were also included:
And build up to the figure above.
timhoffm commentedon May 17, 2024
It is difficult, but better than nothing at all. We don't have the capacity to pamper users for every edge case.
I'm -0.5 on a tutorial. First, I'd like to have both the separate plane and the intersecting plane show up in the gallery overview. Second, going into an extensive tutorial for working around limitations and insufficiencies of our API feels the wrong way. It would be better to invest that time into improving the API, i.e. the quadmesh-like parallel-axes plane.
But as always, if you are convinced that a tutorial is worth the effort, you are free to create one.
story645 commentedon May 17, 2024
Maybe the solution is more code comments and not wrapping the figure code in figure_3D_array_slices, but my concern is that the example as written might be hard to modify and cross apply to a slightly different use case, which I think is the core purpose of examples and not an edge case. And a very consistent complaint about the examples is that they tend to do too much in a way where it's hard to see what code is responsible for which aspect of the visualization.
I think of it more as a tutorial showing how to use the existing API to build something that the API doesn't provide out of the box, but yeah if a new function makes more sense than I'm not opposed to that.