#!/usr/bin/env python3
"""Create a 4096x4096 allRGB image resembling a light micrograph of dopamine crystals.
The generated image contains every possible 24-bit RGB color exactly once
(16,777,216 pixels total), while trying to visually evoke the colourful,
interference-like crystal patterns seen in polarized/differential-interference
light micrographs.
Requires: numpy, pillow, scipy
Example:
python dopamine_allrgb_crystals.py --output dopamine_crystals_allrgb.png --seed 17
Notes:
- The allRGB constraint forces the final image to be 4096 x 4096.
- This is an artistic approximation, not a scientific reconstruction.
- Use --save-scene if you want to inspect the pre-allRGB crystal scene.
"""
from __future__ import annotations
import argparse
import gc
import math
from pathlib import Path
import numpy as np
from PIL import Image, ImageFilter
from scipy.ndimage import gaussian_filter
SIZE = 4096
SCENE_SIZE = 1024
N = SIZE * SIZE
def hsv_to_rgb(h: np.ndarray, s: np.ndarray, v: np.ndarray) -> np.ndarray:
"""Vectorized HSV -> RGB, inputs in [0,1]. Returns uint8 RGB."""
h = np.mod(h, 1.0)
s = np.clip(s, 0.0, 1.0)
v = np.clip(v, 0.0, 1.0)
i = np.floor(h * 6.0).astype(np.int32)
f = h * 6.0 - i
p = v * (1.0 - s)
q = v * (1.0 - f * s)
t = v * (1.0 - (1.0 - f) * s)
i %= 6
r = np.empty_like(v)
g = np.empty_like(v)
b = np.empty_like(v)
m = i == 0
r[m], g[m], b[m] = v[m], t[m], p[m]
m = i == 1
r[m], g[m], b[m] = q[m], v[m], p[m]
m = i == 2
r[m], g[m], b[m] = p[m], v[m], t[m]
m = i == 3
r[m], g[m], b[m] = p[m], q[m], v[m]
m = i == 4
r[m], g[m], b[m] = t[m], p[m], v[m]
m = i == 5
r[m], g[m], b[m] = v[m], p[m], q[m]
rgb = np.stack((r, g, b), axis=-1)
return np.clip(np.round(rgb * 255.0), 0, 255).astype(np.uint8)
def morton_encode_rgb(r: np.ndarray, g: np.ndarray, b: np.ndarray) -> np.ndarray:
"""Encode 8-bit R,G,B arrays into a 24-bit Morton/Z-order key."""
r = r.astype(np.uint32, copy=False)
g = g.astype(np.uint32, copy=False)
b = b.astype(np.uint32, copy=False)
code = np.zeros(r.shape, dtype=np.uint32)
for bit in range(8):
code |= ((r >> bit) & 1) << (3 * bit + 2)
code |= ((g >> bit) & 1) << (3 * bit + 1)
code |= ((b >> bit) & 1) << (3 * bit)
return code
def morton_decode_rgb(code: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Decode 24-bit Morton/Z-order keys back into 8-bit R,G,B arrays."""
code = code.astype(np.uint32, copy=False)
r = np.zeros(code.shape, dtype=np.uint32)
g = np.zeros(code.shape, dtype=np.uint32)
b = np.zeros(code.shape, dtype=np.uint32)
for bit in range(8):
r |= ((code >> (3 * bit + 2)) & 1) << bit
g |= ((code >> (3 * bit + 1)) & 1) << bit
b |= ((code >> (3 * bit + 0)) & 1) << bit
return r.astype(np.uint8), g.astype(np.uint8), b.astype(np.uint8)
def make_crystal_scene(seed: int, small: int = SCENE_SIZE) -> np.ndarray:
rng = np.random.default_rng(seed)
h = w = small
y, x = np.mgrid[0:h, 0:w].astype(np.float32)
# Soft colour fog and uneven illumination as a base.
noise_large = gaussian_filter(rng.normal(size=(h, w)).astype(np.float32), sigma=24)
noise_mid = gaussian_filter(rng.normal(size=(h, w)).astype(np.float32), sigma=8)
illum = 0.55 + 0.18 * noise_large / (np.std(noise_large) + 1e-6)
illum = np.clip(illum, 0.18, 0.95)
# Seed many crystal domains with differing orientations and colour behaviour.
n_domains = 95
centers_x = rng.uniform(-60, w + 60, n_domains).astype(np.float32)
centers_y = rng.uniform(-60, h + 60, n_domains).astype(np.float32)
theta = rng.uniform(0, 2 * np.pi, n_domains).astype(np.float32)
stretch_u = rng.uniform(0.65, 1.8, n_domains).astype(np.float32)
stretch_v = rng.uniform(0.65, 1.8, n_domains).astype(np.float32)
base_hue = rng.uniform(0.0, 1.0, n_domains).astype(np.float32)
fan_count = rng.integers(5, 16, n_domains).astype(np.float32)
band_freq = rng.uniform(0.055, 0.19, n_domains).astype(np.float32)
radial_twist = rng.uniform(0.01, 0.09, n_domains).astype(np.float32)
phase1 = rng.uniform(0, 2 * np.pi, n_domains).astype(np.float32)
phase2 = rng.uniform(0, 2 * np.pi, n_domains).astype(np.float32)
sat_bias = rng.uniform(-0.05, 0.32, n_domains).astype(np.float32)
val_bias = rng.uniform(-0.08, 0.22, n_domains).astype(np.float32)
neutralish = rng.random(n_domains) < 0.28
best = np.full((h, w), np.inf, dtype=np.float32)
second = np.full((h, w), np.inf, dtype=np.float32)
idx = np.zeros((h, w), dtype=np.int16)
print('Building crystal domains...', flush=True)
for i in range(n_domains):
dx = x - centers_x[i]
dy = y - centers_y[i]
c = math.cos(float(theta[i]))
s = math.sin(float(theta[i]))
u = (dx * c + dy * s) / stretch_u[i]
v = (-dx * s + dy * c) / stretch_v[i]
# Add a tiny directional warp so some domains become fan/plate-like.
d = u * u + v * v + 0.0018 * (u * v)
m = d < best
second = np.where(m, best, second)
idx = np.where(m, i, idx)
best = np.where(m, d, best)
m2 = (~m) & (d < second)
second = np.where(m2, d, second)
domain = idx
cx = centers_x[domain]
cy = centers_y[domain]
dx = x - cx
dy = y - cy
radius = np.sqrt(np.maximum(best, 0.0))
angle = np.arctan2(dy, dx)
edge_dist = np.sqrt(np.maximum(second - best, 0.0))
fan = 0.5 + 0.5 * np.sin(
fan_count[domain] * angle
+ radius * band_freq[domain] * 7.5
+ phase1[domain]
+ radial_twist[domain] * radius * angle
)
bands = 0.5 + 0.5 * np.sin(radius * (4.0 * band_freq[domain]) + phase2[domain])
feather = 0.5 + 0.5 * np.sin(
0.7 * fan_count[domain] * angle - radius * band_freq[domain] * 6.5 + phase1[domain] * 0.55
)
shimmer = 0.5 + 0.5 * np.sin(0.022 * x + 0.031 * y + 3.2 * noise_mid)
boundary = np.clip(edge_dist / 5.5, 0.0, 1.0)
chips = gaussian_filter((rng.random((h, w)) > 0.985).astype(np.float32), sigma=1.0)
chips = np.clip(chips * 3.5, 0.0, 1.0)
hue = (
base_hue[domain]
+ 0.18 * (fan - 0.5)
+ 0.12 * (bands - 0.5)
+ 0.06 * (shimmer - 0.5)
+ 0.03 * noise_mid
) % 1.0
sat = (
0.40
+ sat_bias[domain]
+ 0.35 * fan
+ 0.12 * feather
+ 0.08 * shimmer
)
val = (
0.30
+ val_bias[domain]
+ 0.34 * bands
+ 0.18 * feather
+ 0.08 * shimmer
)
# Some domains become more silvery / pastel, like real crystal regions.
sat = np.where(neutralish[domain], sat * 0.42, sat)
val = np.where(neutralish[domain], val + 0.12, val)
sat *= 0.78 + 0.22 * boundary
val *= 0.72 + 0.28 * boundary
val = np.clip(val + 0.18 * chips, 0.0, 1.0)
# Add a few stronger plate-like highlights and darker creases.
sparkle = gaussian_filter(rng.random((h, w)).astype(np.float32), sigma=2.2)
sparkle = np.clip((sparkle - 0.72) * 2.6, 0.0, 1.0)
crease = np.clip(1.0 - boundary, 0.0, 1.0)
sat = np.clip(sat, 0.04, 1.0)
val = np.clip(val * illum + 0.10 * sparkle - 0.12 * crease, 0.0, 1.0)
rgb = hsv_to_rgb(hue, sat, val)
# Gentle blend to mimic microscope optics while keeping crystal edges.
img = Image.fromarray(rgb, 'RGB').filter(ImageFilter.GaussianBlur(radius=0.45))
img = img.resize((SIZE, SIZE), Image.Resampling.LANCZOS)
# Final colour/grain treatment at full size.
arr = np.asarray(img, dtype=np.float32)
fine = rng.normal(0.0, 5.0, size=arr.shape).astype(np.float32)
broad = gaussian_filter(rng.normal(size=(SIZE, SIZE)).astype(np.float32), sigma=18)
broad = broad[..., None]
arr += fine + 7.5 * broad / (np.std(broad) + 1e-6)
# Mild channel cross-variation for iridescent interference feel.
arr[..., 0] *= 1.03
arr[..., 1] *= 0.995
arr[..., 2] *= 1.02
arr = np.clip(arr, 0, 255).astype(np.uint8)
return arr
def save_scene_preview(scene: np.ndarray, path: str | Path) -> None:
Image.fromarray(scene, 'RGB').save(path)
print(f'Saved scene preview: {path}', flush=True)
def render_allrgb(seed: int, output: str | Path, save_scene: str | None = None, chunk: int = 500_000) -> None:
scene = make_crystal_scene(seed)
if save_scene:
save_scene_preview(scene, save_scene)
print('Computing colour-order keys for target scene...', flush=True)
r = scene[..., 0].ravel()
g = scene[..., 1].ravel()
b = scene[..., 2].ravel()
keys = morton_encode_rgb(r, g, b)
del scene, r, g, b
gc.collect()
print('Sorting pixels by target colour structure...', flush=True)
pixel_order = np.argsort(keys, kind='stable')
del keys
gc.collect()
print('Assigning every RGB value exactly once...', flush=True)
out = np.empty((N, 3), dtype=np.uint8)
for start in range(0, N, chunk):
end = min(start + chunk, N)
codes = np.arange(start, end, dtype=np.uint32)
rr, gg, bb = morton_decode_rgb(codes)
out[pixel_order[start:end], 0] = rr
out[pixel_order[start:end], 1] = gg
out[pixel_order[start:end], 2] = bb
if start == 0 or end == N or (start // chunk) % 8 == 0:
print(f' {end:,} / {N:,} colours assigned', flush=True)
del pixel_order
gc.collect()
Image.fromarray(out.reshape(SIZE, SIZE, 3), 'RGB').save(output, optimize=False)
print(f'Saved {output}: {SIZE} x {SIZE}, all {N:,} RGB colours exactly once.', flush=True)
def main() -> None:
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--output', default='dopamine_crystals_allrgb.png', help='Final allRGB PNG output')
parser.add_argument('--seed', type=int, default=17, help='Random seed for composition')
parser.add_argument('--save-scene', default=None, help='Optional path to save the pre-allRGB crystal scene')
parser.add_argument('--chunk', type=int, default=500_000, help='Chunk size for allRGB assignment')
args = parser.parse_args()
render_allrgb(seed=args.seed, output=args.output, save_scene=args.save_scene, chunk=args.chunk)
if __name__ == '__main__':
main()
