3D slice of a 4D surface with Python

slice_Tiger_01.png

slice_Tiger_01.py

# -*- 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)

slice_Tiger_01_anim.py

# -*- 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

slicedTiger.gif