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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
"""
.. _marching_cubes_example:
Marching Cubes
~~~~~~~~~~~~~~
Generate a surface from a scalar field using the flying edges and
marching cubes filters as provided by the :func:`contour
<pyvista.DataSetFilters.contour>` filter.
Special thanks to GitHub user `stla <https://gist.github.com/stla>`_
for providing examples.
"""
from __future__ import annotations
import numpy as np
import pyvista as pv
# %%
# Spider Cage
# ~~~~~~~~~~~
# Use the marching cubes algorithm to extract the isosurface
# generated from the spider cage function.
a = 0.9
def spider_cage(x, y, z):
x2 = x * x
y2 = y * y
x2_y2 = x2 + y2
return (np.sqrt((x2 - y2) ** 2 / x2_y2 + 3 * (z * np.sin(a)) ** 2) - 3) ** 2 + 6 * (
np.sqrt((x * y) ** 2 / x2_y2 + (z * np.cos(a)) ** 2) - 1.5
) ** 2
# create a uniform grid to sample the function with
n = 100
x_min, y_min, z_min = -5, -5, -3
grid = pv.ImageData(
dimensions=(n, n, n),
spacing=(abs(x_min) / n * 2, abs(y_min) / n * 2, abs(z_min) / n * 2),
origin=(x_min, y_min, z_min),
)
x, y, z = grid.points.T
# sample and plot
values = spider_cage(x, y, z)
mesh = grid.contour([1], values, method='marching_cubes')
dist = np.linalg.norm(mesh.points, axis=1)
mesh.plot(scalars=dist, smooth_shading=True, cmap='plasma', show_scalar_bar=False)
# %%
# Barth Sextic
# ~~~~~~~~~~~~
# Use the flying edges algorithm to extract the isosurface
# generated from the Barth sextic function.
phi = (1 + np.sqrt(5)) / 2
phi2 = phi * phi
def barth_sextic(x, y, z):
x2 = x * x
y2 = y * y
z2 = z * z
arr = (
3 * (phi2 * x2 - y2) * (phi2 * y2 - z2) * (phi2 * z2 - x2)
- (1 + 2 * phi) * (x2 + y2 + z2 - 1) ** 2
)
nan_mask = x2 + y2 + z2 > 3.1
arr[nan_mask] = np.nan
return arr
# create a uniform grid to sample the function with
n = 100
k = 2.0
x_min, y_min, z_min = -k, -k, -k
grid = pv.ImageData(
dimensions=(n, n, n),
spacing=(abs(x_min) / n * 2, abs(y_min) / n * 2, abs(z_min) / n * 2),
origin=(x_min, y_min, z_min),
)
x, y, z = grid.points.T
# sample and plot
values = barth_sextic(x, y, z)
mesh = grid.contour([0], values, method='flying_edges')
dist = np.linalg.norm(mesh.points, axis=1)
mesh.plot(scalars=dist, smooth_shading=True, cmap='plasma', show_scalar_bar=False)
# %%
# Animate Barth Sextic
# ~~~~~~~~~~~~~~~~~~~~
# Show 20 frames of various isocurves extracted from the Barth sextic
# function.
def angle_to_range(angle):
return -2 * np.sin(angle)
pl = pv.Plotter(window_size=[800, 800], off_screen=True)
pl.open_gif('barth_sextic.gif')
for angle in np.linspace(0, np.pi, 20, endpoint=False):
# clear the plotter before adding each frame's mesh
pl.clear()
pl.enable_lightkit()
mesh = grid.contour([angle_to_range(angle)], values, method='flying_edges')
dist = np.linalg.norm(mesh.points, axis=1)
pl.add_mesh(
mesh,
scalars=dist,
smooth_shading=True,
rng=[0.5, 1.5],
cmap='plasma',
show_scalar_bar=False,
)
pl.write_frame()
pl.close()
# %%
# .. tags:: filter
|
