Sierpinski Pyramid Voxels

import matplotlib.pyplot as plt
import numpy as np

def sierpinski_3d(level, size):
    if level == 0:
        return np.ones((size, size, size), dtype=bool)
    
    half = size // 2
    sub = sierpinski_3d(level - 1, half)
    result = np.zeros((size, size, size), dtype=bool)
    
    # 4 corner tetrahedra
    result[:half, :half, :half] = sub
    result[half:, :half, :half] = sub
    result[:half, half:, :half] = sub
    result[half:size//2+half//2, half:size//2+half//2, half:] = sub[:half//2, :half//2, :half]
    
    return result

n = 16
voxels = np.zeros((n, n, n), dtype=bool)

# Simple pyramid approximation
for z in range(n):
    margin = z // 2
    if margin < n // 2:
        voxels[margin:n-margin, margin:n-margin, z] = True

# Gradient from dark_red #9C2007 to red #C82909 to orange #F54927
z_idx = np.indices((n, n, n))[2]
colors = np.empty(voxels.shape + (4,), dtype=np.float32)
colors[..., 3] = 0

norm_z = z_idx / (n - 1)
colors[voxels, 0] = 0.61 + 0.35 * norm_z[voxels]  # 9C to F5
colors[voxels, 1] = 0.13 + 0.16 * norm_z[voxels]  # 20 to 49
colors[voxels, 2] = 0.03 + 0.12 * norm_z[voxels]  # 07 to 27
colors[voxels, 3] = 0.85

fig = plt.figure(figsize=(10, 8), facecolor='#020B14')
ax = fig.add_subplot(111, projection='3d', facecolor='#020B14')

ax.voxels(voxels, facecolors=colors, edgecolors='#1e293b', linewidth=0.2)

ax.set_xlabel('X', fontsize=11, color='#94a3b8', labelpad=10)
ax.set_ylabel('Y', fontsize=11, color='#94a3b8', labelpad=10)
ax.set_zlabel('Z', fontsize=11, color='#94a3b8', labelpad=10)
ax.set_title('Sierpinski Pyramid Voxels', fontsize=14, color='white', fontweight='bold', pad=20)

ax.tick_params(colors='#64748b', labelsize=8)
ax.xaxis.pane.fill = False
ax.yaxis.pane.fill = False
ax.zaxis.pane.fill = False
ax.xaxis.pane.set_edgecolor('#1e293b')
ax.yaxis.pane.set_edgecolor('#1e293b')
ax.zaxis.pane.set_edgecolor('#1e293b')
ax.grid(True, alpha=0.2, color='#475569')

plt.tight_layout()
plt.show()