1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
"""
===================
Intersecting planes
===================
This examples demonstrates drawing intersecting planes in 3D. It is a generalization
of :doc:`/gallery/mplot3d/imshow3d`.
Drawing intersecting planes in `.mplot3d` is complicated, because `.mplot3d` is not a
real 3D renderer, but only projects the Artists into 3D and draws them in the right
order. This does not work correctly if Artists overlap each other mutually. In this
example, we lift the problem of mutual overlap by segmenting the planes at their
intersections, making four parts out of each plane.
This examples only works correctly for planes that cut each other in haves. This
limitation is intentional to keep the code more readable. Cutting at arbitrary
positions would of course be possible but makes the code even more complex.
Thus, this example is more a demonstration of the concept how to work around
limitations of the 3D visualization, it's not a refined solution for drawing
arbitrary intersecting planes, which you can copy-and-paste as is.
"""
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 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()
# %%
# .. tags::
# plot-type: 3D,
# level: advanced
|
