Skip to content

Instantly share code, notes, and snippets.

@stla
Last active Oct 18, 2021
Embed
What would you like to do?
3D slice of a 4D surface with Python
# -*- coding: utf-8 -*-
import numpy as np
import pyvista as pv
R1 = 2
R2 = 2
r = 0.5
def s(u, v, w):
return np.array([
np.cos(u) * (R1 + r*np.cos(w)),
np.sin(u) * (R1 + r*np.cos(w)),
np.cos(v) * (R2 + r*np.sin(w)),
np.sin(v) * (R2 + r*np.sin(w))
])
a = np.array([1, 1, 1, 1])
b = 2 # plane x+y+z+w = 2
x0 = np.array([b, b, b, b]) / 4 # a point in this plane
nrml = a / np.linalg.norm(a) # unit normal
def f(u, v, w):
return np.matmul(np.transpose(s(u, v, w)), nrml.reshape((4,1)))
u, v, w = np.mgrid[0:2*np.pi:100j, 0:2*np.pi:100j, 0:2*np.pi:100j]
values = f(np.transpose(u), np.transpose(v), np.transpose(w))
grid = pv.StructuredGrid(u, v, w)
grid.point_data["values"] = values.ravel(order="F")
isosurf = grid.contour(isosurfaces=[np.vdot(x0, nrml)])
mesh0 = isosurf.extract_geometry()
vs0 = mesh0.points
vs = s(vs0[:,0], vs0[:,1], vs0[:,2])
def rotationMatrix4D(v1, v2):
v1 = v1 / np.linalg.norm(v1)
v2 = v2 / np.linalg.norm(v2)
v = (v1 + v2).reshape((4,1))
return 2*np.matmul(v, np.transpose(v)) / (np.vdot(v,v)) - np.eye(4)
Rot = rotationMatrix4D(nrml, np.array([0,0,0,1]))
VSprime = np.matmul(Rot, vs)[[0,1,2], :]
mesh = pv.PolyData(np.transpose(VSprime), mesh0.faces)
mesh.plot(smooth_shading=True, color="orange", specular=5)
# -*- coding: utf-8 -*-
import numpy as np
import pyvista as pv
R1 = 3
R2 = 2
r = 0.25
def s(u, v, w):
return np.array([
np.cos(u) * (R1 + r*np.cos(w)),
np.sin(u) * (R1 + r*np.cos(w)),
np.cos(v) * (R2 + r*np.sin(w)),
np.sin(v) * (R2 + r*np.sin(w))
])
a = np.array([1, 1, 1, 1]) # plane x+y+z+w = b
nrml = a / np.linalg.norm(a) # unit normal
def f(u, v, w):
return np.matmul(np.transpose(s(u, v, w)), nrml.reshape((4,1)))
u, v, w = np.mgrid[0:2*np.pi:150j, 0:2*np.pi:150j, 0:2*np.pi:150j]
values = f(np.transpose(u), np.transpose(v), np.transpose(w))
grid = pv.StructuredGrid(u, v, w)
grid.point_data["values"] = values.ravel(order="F")
def rotationMatrix4D(v1, v2):
v1 = v1 / np.linalg.norm(v1)
v2 = v2 / np.linalg.norm(v2)
v = (v1 + v2).reshape((4,1))
return 2*np.matmul(v, np.transpose(v)) / (np.vdot(v,v)) - np.eye(4)
Rot = rotationMatrix4D(nrml, np.array([0,0,0,1]))
b_ = np.linspace(-5, 5, 60)
for i, b in enumerate(b_):
print("i", i)
i = str(i+1) if i >= 9 else "0" + str(i+1)
pngfile = "zpics/tiger" + i + ".png"
x0 = np.array([b, b, b, b]) / 4 # a point in the plane
isosurf = grid.contour(isosurfaces=[np.vdot(x0, nrml)])
mesh0 = isosurf.extract_geometry()
vs0 = mesh0.points
vs = s(vs0[:,0], vs0[:,1], vs0[:,2])
VSprime = np.matmul(Rot, vs)[[0,1,2], :]
mesh = pv.PolyData(np.transpose(VSprime), mesh0.faces)
pl = pv.Plotter(off_screen=True, window_size=[500, 500])
pl.set_background("#363940")
pl.add_mesh(mesh, smooth_shading=True, color="orange", specular=5)
pl.show(screenshot=pngfile)
# magick convert -dispose previous -loop 0 -delay 10 -duplicate 1,-2-1 tiger*.png slicedTiger.gif
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment