import numpy as np
import colorsys
import math
from itertools import product
from PIL import Image
# === Configuration ===
SIZE = 4096
RADIUS = 6 # dot radius in pixels
ALPHA_BG = 128 # 50% opacity for background
ALPHA_FG = 255 # fully opaque for lines & dots
TOTAL_PIXELS = SIZE * SIZE
TOTAL_COLORS = 2**24 # 16,777,216
# === 1) Generate the 240 E₈ roots in ℝ⁸ ===
roots = []
# A-type: ±e_i ± e_j
for i in range(8):
for j in range(i+1, 8):
v = np.zeros(8, dtype=float)
v[i], v[j] = 1, -1
roots.append(v.copy()); roots.append((-v).copy())
v[i], v[j] = 1, 1
roots.append(v.copy()); roots.append((-v).copy())
# B-type: (±½,…,±½) with product(signs) = +1
for signs in product([1, -1], repeat=8):
if np.prod(signs) == 1:
roots.append(np.array(signs, dtype=float) * 0.5)
roots = np.vstack(roots) # shape (240,8)
# === 2) Project into 2D via PCA/SVD ===
roots_centered = roots - roots.mean(axis=0)
U, S, Vh = np.linalg.svd(roots_centered, full_matrices=False)
proj2d = roots_centered @ Vh.T[:, :2] # (240,2)
# === 3) Scale & translate into pixel coordinates ===
marg = 0.05 * SIZE
scale = (SIZE/2 - marg) / np.max(np.linalg.norm(proj2d, axis=1))
coords_f = proj2d * scale + SIZE/2 # floats
coords = np.round(coords_f).astype(int)
# === 4) Find all edges (⟨α,β⟩ = 1) ===
gram = roots @ roots.T
i_u, j_u = np.where(np.triu(np.isclose(gram, 1.0), k=1))
edges = list(zip(i_u.tolist(), j_u.tolist()))
# === 5) Bresenham line algorithm ===
def bresenham(x0, y0, x1, y1):
dx = abs(x1 - x0); sx = 1 if x0 < x1 else -1
dy = abs(y1 - y0); sy = 1 if y0 < y1 else -1
err = dx - dy
pixels = []
while True:
pixels.append((x0, y0))
if x0 == x1 and y0 == y1:
break
e2 = 2 * err
if e2 > -dy:
err -= dy
x0 += sx
if e2 < dx:
err += dx
y0 += sy
return pixels
# === 6) Generate circle mask offsets for dots ===
circle_offsets = []
for dy in range(-RADIUS, RADIUS + 1):
dx_max = int(math.floor(math.sqrt(RADIUS*RADIUS - dy*dy)))
for dx in range(-dx_max, dx_max + 1):
circle_offsets.append((dx, dy))
# === 7) Collect pixel coordinates for lines & dots ===
used = np.zeros((SIZE, SIZE), dtype=bool)
line_pixels = []
for i, j in edges:
x0, y0 = coords[i]; x1, y1 = coords[j]
for (x, y) in bresenham(x0, y0, x1, y1):
if not used[y, x]:
line_pixels.append((x, y))
used[y, x] = True
dot_pixels = []
for x0, y0 in coords:
for dx, dy in circle_offsets:
x, y = x0 + dx, y0 + dy
if 0 <= x < SIZE and 0 <= y < SIZE and not used[y, x]:
dot_pixels.append((x, y))
used[y, x] = True
N_line = len(line_pixels)
N_dot = len(dot_pixels)
N_bg = TOTAL_PIXELS - N_line - N_dot
assert N_bg + N_line + N_dot == TOTAL_PIXELS, "Pixel count mismatch!"
# === 8) Build & sort full 24-bit palette by HSV (v, s descending) ===
# Create array of shape (16_777_216, 3)
palette = np.empty((TOTAL_COLORS, 3), dtype=np.uint8)
idxs = np.arange(TOTAL_COLORS, dtype=np.uint32)
palette[:, 0] = (idxs >> 16) & 0xFF
palette[:, 1] = (idxs >> 8) & 0xFF
palette[:, 2] = (idxs ) & 0xFF
# Convert to HSV
# scaled to [0,1] floats
rgb_norm = palette.astype(float) / 255.0
hsv = np.array([colorsys.rgb_to_hsv(*c) for c in rgb_norm])
v = hsv[:, 2]
s = hsv[:, 1]
# Sort so brightest & most saturated first:
order = np.lexsort((-(s), -(v)))
sorted_palette = palette[order]
# === 9) Assign colors to lines, dots, and background ===
line_colors = sorted_palette[:N_line]
dot_colors = sorted_palette[N_line:N_line+N_dot]
bg_colors = sorted_palette[N_line+N_dot:]
# === 10) Create RGBA image and paint ===
img = Image.new("RGBA", (SIZE, SIZE))
pixels = img.load()
# Paint lines
for (x, y), (r, g, b) in zip(line_pixels, line_colors):
pixels[x, y] = (int(r), int(g), int(b), ALPHA_FG)
# Paint dots
for (x, y), (r, g, b) in zip(dot_pixels, dot_colors):
pixels[x, y] = (int(r), int(g), int(b), ALPHA_FG)
# Paint background
bg_iter = iter(bg_colors)
for y in range(SIZE):
for x in range(SIZE):
if not used[y, x]:
r, g, b = next(bg_iter)
pixels[x, y] = (int(r), int(g), int(b), ALPHA_BG)
# === 11) Save PNG ===
img.save("E8_Petrie_exact24bit_4096.png")
print("✅ Generated 4096×4096 E₈ Petrie projection using every 24-bit RGB color exactly once.")
