How to use pythreejs to plot a superellipsoid

A superellipsoid is given by a parametric function and the equation is very similar to an ellipse equation. We only have different exponents which give us different shapes. For more informations: https://en.wikipedia.org/wiki/Superellipsoid.

The idea of this example is to construct the mesh of the square $[0, 1]\times[0,1]$ and to do a projection of these points on the superillipse which is the 2D shape and then to do a spherical product to have the 3D shape.


In [2]:
import numpy as np

n = 10 # number of discretisation points for the square in each direction 
x_box = np.concatenate((np.linspace(-1, 1., n), np.ones(n-2), np.linspace(1, -1., n), -np.ones(n-2)))
y_box = np.concatenate((-np.ones(n-1), np.linspace(-1, 1., n), np.ones(n-2), np.linspace(1, -1., n-1, endpoint=False)))
nx_box = x_box.size

coords = np.empty((nx_box**2, 3))

def superellipse(rx, ry, m):
    """
    superellipse formula with the projection of the unit square
    
    Parameters
    ----------
    rx : the radius in the x direction 
    ry : the radius in the y direction 
    m : the exponent of the superellipse
    
    Output
    ------
    the coordinates of the superellipse
    """
    return x_box*rx*(1. - .5*np.abs(y_box)**(2./m))**(m/2.),  y_box*ry*(1. - .5*np.abs(x_box)**(2./m))**(m/2.)

In [3]:
def superellipsoid(rx, ry, rz, m1, m2):
    """
    superellipsoid formula with the spherical product of two superellipse
    and update of the global coords array
    
    Parameters
    ----------
    rx : the radius in the x direction 
    ry : the radius in the y direction 
    rz : the radius in the z direction 
    m1 : the exponent of the first superellipse
    m2 : the exponent of the second superellipse
    """    
    gx, gy = superellipse(1, 1, m2)
    hx, hy = superellipse(1, 1, m1)

    
    coords[:, 0] = rx*(gx[np.newaxis, :]*hx[:, np.newaxis]).flatten()
    coords[:, 1] = ry*(gx[np.newaxis, :]*hy[:, np.newaxis]).flatten()
    coords[:, 2] = rz*(gy[np.newaxis, :]*np.ones(hx.size)[:, np.newaxis]).flatten()

In [4]:
# superellipsoid parameters
rx = ry = rz = 1.
m1 = m2 = 1.

superellipsoid(rx, ry, rz, m1, m2)

We construct the triangulation by using the ConveHull function in scipy.


In [5]:
import scipy.spatial as spatial

cvx = spatial.ConvexHull(coords)

In [6]:
from pythreejs import *
from IPython.display import display

surf_g = PlainGeometry(vertices=coords.tolist(), faces=cvx.simplices.tolist())
surf = Mesh(geometry=surf_g, material=BasicMaterial(color='green', wireframe=True))
scene = Scene(children=[surf, AmbientLight(color='#777777')])
c = PerspectiveCamera(position=[2, 2, 3], up=[0, 0, 1],
                      children=[DirectionalLight(color='white',
                                                 position=[3, 5, 1],
                                                 intensity=0.6)])
renderer = Renderer(camera=c, scene=scene, controls=[OrbitControls(controlling=c)])
display(renderer)

In [7]:
from ipywidgets import FloatSlider, HBox, VBox

m1_slider, m2_slider = (FloatSlider(description='m1', min=0.01, max=4.0, step=0.01, value=m1,
                                            continuous_update=False, orientation='vertical'),
                        FloatSlider(description='m2', min=0.01, max=4.0, step=0.01, value=m2,
                                            continuous_update=False, orientation='vertical'))

In [8]:
rx_slider, ry_slider, rz_slider = (FloatSlider(description='rx', min=0.01, max=10.0, step=0.01, value=rx, 
                                               continuous_update=False, orientation='horizontal'),
                                   FloatSlider(description='ry', min=0.01, max=10.0, step=0.01, value=ry, 
                                               continuous_update=False, orientation='horizontal'),
                                   FloatSlider(description='rz', min=0.01, max=10.0, step=0.01, value=rz, 
                                               continuous_update=False, orientation='horizontal'))

In [9]:
def update(change):
    superellipsoid(rx_slider.value, ry_slider.value, rz_slider.value, 
                   m1_slider.value, m2_slider.value)
    surf_g.vertices = coords.tolist()
    
m1_slider.observe(update, names=['value'])
m2_slider.observe(update, names=['value'])
rx_slider.observe(update, names=['value'])
ry_slider.observe(update, names=['value'])
rz_slider.observe(update, names=['value'])

In [10]:
VBox([HBox([renderer, m1_slider, m2_slider]), rx_slider, ry_slider, rz_slider])

In [ ]: