Pumpkin surface.

pumpkin.png

pumpkin.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Based on Matlab's code

https://www.instagram.com/p/CVkry4ggsmw/

@author: Nicolás Guarín-Zapata
@date: October 2021
"""
import numpy as np
import pyvista as pv


def sphere(n):
    """Generates a unit sphere with a n by n angular grid"""
    phi, theta = np.mgrid[-np.pi:np.pi:n*1j, -np.pi/2:np.pi/2:n*1j]
    x = np.cos(phi) * np.cos(theta)
    y = np.cos(phi) * np.sin(theta)
    z = np.sin(phi)
    return x, y, z


#%% Pumpkin
npts = 500
bumps = 6
bdepth = 0.1
bdepth2 = 0.02
dimple = 0.2
width_r = 1
height_r = 0.8

Xs, Ys, Zs = sphere(npts)
Rxy = -bdepth*(1 - np.mod(np.linspace(0, 2*bumps, npts), 2))**2 \
      -bdepth2 *(1 - np.mod(np.linspace(0, 4*bumps, npts), 2))**2
Rz = dimple*(0 - np.linspace(1, -1, npts)**4)
Rxy, Rz =  np.meshgrid(Rxy, Rz)
Xp = (width_r + Rxy) * Xs
Yp = (width_r + Rxy) * Ys
Zp = (height_r + Rz) * Zs * (Rxy + 1)
Cp = np.hypot(np.hypot(Xp, Yp), width_r * Zs * (Rxy + 1))


#%% Stem
sheight = 0.5
scurve = 0.4
srad = np.outer(np.array([0.1, 0.06]*bumps + [0.1]),
                np.array([1.5, 1] + [0.7]*6))
theta, phi = np.meshgrid(np.linspace(0, np.pi/2, srad.shape[1]),
                         np.linspace(0, 2*np.pi, srad.shape[0]))
Xs = (scurve - np.cos(phi) * srad) * np.cos(theta) - scurve
Zs = (sheight - np.cos(phi) * srad) * np.sin(theta)\
    + height_r - max(0, dimple*0.9)
Ys = -np.sin(phi)* srad


#%% Visualization
pumpkin = pv.StructuredGrid(Xp, Yp, Zp)
stem = pv.StructuredGrid(Xs, Ys, Zs)

plotter = pv.Plotter(window_size=[768, 768], polygon_smoothing=True)
plotter.add_mesh(pumpkin, color="orange", smooth_shading=True,
                 ambient=0.2, diffuse=0.9, specular=0.6, specular_power=80)
plotter.add_mesh(stem, color="green", smooth_shading=True)
plotter.set_background("#333333")
# plotter.add_text("@nicoguaro", position="lower_right")
plotter.show(screenshot="pumpkin.png")

pumpkin_mpl.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Based on Matlab's code

https://www.instagram.com/p/CVkry4ggsmw/

@author: Nicolás Guarín-Zapata
@date: October 2021
"""
import numpy as np
import matplotlib.pyplot as plt


def sphere(n):
    """Generates a unit sphere with a n by n angular grid"""
    phi, theta = np.mgrid[-np.pi:np.pi:n*1j, -np.pi/2:np.pi/2:n*1j]
    x = np.cos(phi) * np.cos(theta)
    y = np.cos(phi) * np.sin(theta)
    z = np.sin(phi)
    return x, y, z


#%% Pumpkin
npts = 200
bumps = 10
bdepth = 0.1
bdepth2 = 0.02
dimple = 0.2
width_r = 1
height_r = 0.8

Xs, Ys, Zs = sphere(npts)
Rxy = -bdepth*(1 - np.mod(np.linspace(0, 2*bumps, npts), 2))**2 \
      -bdepth2 *(1 - np.mod(np.linspace(0, 4*bumps, npts), 2))**2
Rz = dimple*(0 - np.linspace(1, -1, npts)**4)
Rxy, Rz =  np.meshgrid(Rxy, Rz)
Xp = (width_r + Rxy) * Xs
Yp = (width_r + Rxy) * Ys
Zp = (height_r + Rz) * Zs * (Rxy + 1)
Cp = np.hypot(np.hypot(Xp, Yp), width_r * Zs * (Rxy + 1))


#%% Stem
sheight = 0.5
scurve = 0.4
srad = np.outer(np.array([0.1, 0.06]*bumps + [0.1]),
                np.array([1.5, 1] + [0.7]*6))
theta, phi = np.meshgrid(np.linspace(0, np.pi/2, srad.shape[1]),
                         np.linspace(0, 2*np.pi, srad.shape[0]))
Xs = (scurve - np.cos(phi) * srad) * np.cos(theta) - scurve
Zs = (sheight - np.cos(phi) * srad) * np.sin(theta)\
    + height_r - max(0, dimple*0.9)
Ys = -np.sin(phi)* srad


#%% Visualization
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
stem = ax.plot_surface(Xs.T, Ys.T, Zs.T, rstride=1, cstride=1, color="green")
pumpkin = ax.plot_surface(Xp, Yp, Zp, rstride=1, cstride=1, color="orange")