#!/usr/bin/env python3
"""Create a 4096x4096 allRGB image resembling a pastel light micrograph of dopamine crystals.
This version is tuned to be closer to the provided reference image:
- softer pastel / pearly palette
- more faceted, plate-like crystal regions
- less swirly than the first crystal version
- still a strict allRGB image: every 24-bit RGB colour appears exactly once
Requires: numpy, pillow, scipy
Example:
python dopamine_allrgb_crystals_v2.py --output dopamine_crystals_allrgb_v2.png --seed 17
Optional:
python dopamine_allrgb_crystals_v2.py --save-scene scene_preview.png
"""
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 = 1152
N = SIZE * SIZE
def hsv_to_rgb(h: np.ndarray, s: np.ndarray, v: np.ndarray) -> np.ndarray:
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]
return np.clip(np.round(np.stack((r, g, b), axis=-1) * 255.0), 0, 255).astype(np.uint8)
def morton_encode_rgb(r: np.ndarray, g: np.ndarray, b: np.ndarray) -> np.ndarray:
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]:
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 polygon_mask_bbox(height: int, width: int, points: np.ndarray):
"""Rasterize polygon into a bounding-box-sized boolean mask."""
min_x = max(0, int(np.floor(points[:, 0].min())) - 2)
max_x = min(width, int(np.ceil(points[:, 0].max())) + 3)
min_y = max(0, int(np.floor(points[:, 1].min())) - 2)
max_y = min(height, int(np.ceil(points[:, 1].max())) + 3)
if min_x >= max_x or min_y >= max_y:
return None
from PIL import ImageDraw
mask_img = Image.new('L', (max_x - min_x, max_y - min_y), 0)
draw = ImageDraw.Draw(mask_img)
pts = [(float(x - min_x), float(y - min_y)) for x, y in points]
draw.polygon(pts, fill=255)
mask = np.asarray(mask_img, dtype=np.uint8) > 0
return min_x, max_x, min_y, max_y, mask
def make_crystal_scene(seed: int, small: int = SCENE_SIZE) -> np.ndarray:
rng = np.random.default_rng(seed)
h = w = small
yy, xx = np.mgrid[0:h, 0:w].astype(np.float32)
# Light microscope paper-white base with slight shading.
base = np.ones((h, w, 3), dtype=np.float32) * 238.0
large_noise = gaussian_filter(rng.normal(size=(h, w)).astype(np.float32), sigma=24)
illum = 6.0 * large_noise / (np.std(large_noise) + 1e-6)
for c in range(3):
base[..., c] += illum
# Build many overlapping faceted crystal plates.
n_crystals = 180
coverage = np.zeros((h, w), dtype=np.float32)
pastel_hues = np.array([
0.00, 0.04, 0.08, 0.12, 0.16, 0.22, 0.30, 0.44, 0.52, 0.58, 0.67, 0.75, 0.83, 0.92
], dtype=np.float32)
print('Building faceted crystal scene...', flush=True)
for _ in range(n_crystals):
cx = rng.uniform(-60, w + 60)
cy = rng.uniform(-60, h + 60)
sides = int(rng.integers(4, 8))
radius = rng.uniform(45, 180)
angle0 = rng.uniform(0, 2 * np.pi)
# Irregular polygon for crystal plate.
angles = np.sort(rng.uniform(0, 2 * np.pi, sides)) + angle0
radii = radius * rng.uniform(0.65, 1.25, sides)
pts = np.stack([cx + np.cos(angles) * radii, cy + np.sin(angles) * radii], axis=1).astype(np.float32)
bbox_data = polygon_mask_bbox(h, w, pts)
if bbox_data is None:
continue
min_x, max_x, min_y, max_y, mask = bbox_data
if mask.sum() < 150:
continue
# Local coordinates and crystal orientation.
theta = rng.uniform(0, 2 * np.pi)
c = math.cos(theta)
s = math.sin(theta)
xloc = xx[min_y:max_y, min_x:max_x]
yloc = yy[min_y:max_y, min_x:max_x]
dx = xloc - cx
dy = yloc - cy
u = dx * c + dy * s
v = -dx * s + dy * c
# Distance to center and to edges.
rr = np.sqrt(dx * dx + dy * dy)
# Use blur of mask to approximate plate depth and edge glow.
mask_f = mask.astype(np.float32)
soft = gaussian_filter(mask_f, sigma=rng.uniform(2.0, 5.0))
edge = np.clip((soft - gaussian_filter(mask_f, sigma=8.0)) * 8.0, -1.0, 1.0)
edge = np.abs(edge)
# Faceted sub-planes: combine 2-3 directional band systems.
f1 = rng.uniform(0.035, 0.095)
f2 = rng.uniform(0.045, 0.13)
f3 = rng.uniform(0.025, 0.085)
phi1 = rng.uniform(0, 2 * np.pi)
phi2 = rng.uniform(0, 2 * np.pi)
phi3 = rng.uniform(0, 2 * np.pi)
bands1 = 0.5 + 0.5 * np.sin(u * f1 * 10.0 + phi1)
bands2 = 0.5 + 0.5 * np.sin((0.72 * u + 0.95 * v) * f2 * 10.0 + phi2)
bands3 = 0.5 + 0.5 * np.sin((1.15 * v - 0.38 * u) * f3 * 10.0 + phi3)
facets = 0.50 * bands1 + 0.32 * bands2 + 0.18 * bands3
facets = gaussian_filter(facets.astype(np.float32), sigma=0.75)
# Radial bursts like interference fans, but much softer than v1.
ang = np.arctan2(dy, dx)
burst_count = rng.integers(6, 18)
burst = 0.5 + 0.5 * np.sin(burst_count * ang + rr * rng.uniform(0.06, 0.16) + rng.uniform(0, 2*np.pi))
burst *= 0.45
# Colour strategy: pastel pearly regions with local rainbow shifts.
base_h = rng.choice(pastel_hues)
hue = (base_h + 0.10 * (facets - 0.5) + 0.08 * (burst - 0.22)) % 1.0
sat = 0.18 + rng.uniform(0.06, 0.22) + 0.26 * facets + 0.12 * burst
sat *= rng.uniform(0.75, 1.10)
sat = np.clip(sat, 0.06, 0.72)
val = 0.70 + rng.uniform(0.02, 0.16) + 0.18 * facets + 0.14 * burst
val += 0.12 * np.clip(soft, 0.0, 1.0)
val -= 0.10 * edge
val = np.clip(val, 0.22, 1.0)
rgb = hsv_to_rgb(hue, sat, val).astype(np.float32)
# Blend onto base with semi-transparency so plates stack like the reference.
alpha = np.clip(0.34 + 0.42 * soft + 0.18 * facets, 0.0, 0.92)
alpha *= mask_f
alpha *= rng.uniform(0.68, 0.92)
sub = base[min_y:max_y, min_x:max_x]
for cch in range(3):
sub[..., cch] = np.where(
mask,
sub[..., cch] * (1.0 - alpha) + rgb[..., cch] * alpha,
sub[..., cch],
)
base[min_y:max_y, min_x:max_x] = sub
coverage[min_y:max_y, min_x:max_x] += alpha
# Add occasional fine crystal fractures/highlights.
if rng.random() < 0.55:
ridge = np.abs(np.sin((u * rng.uniform(0.08, 0.20) + v * rng.uniform(0.02, 0.08) + rng.uniform(0, 2*np.pi))))
ridge = np.clip((ridge - 0.93) * 14.0, 0.0, 1.0) * mask_f
ridge = gaussian_filter(ridge.astype(np.float32), sigma=0.55)
sub = base[min_y:max_y, min_x:max_x]
for cch in range(3):
sub[..., cch] += ridge * rng.uniform(12, 28)
base[min_y:max_y, min_x:max_x] = sub
# Background should remain pale where crystals do not cover.
uncovered = np.clip(1.0 - coverage / max(coverage.max(), 1.0), 0.0, 1.0)
base[..., 0] += uncovered * 10.0
base[..., 1] += uncovered * 8.0
base[..., 2] += uncovered * 11.0
# Gentle microscope diffusion and grain.
base = np.clip(base, 0, 255)
img = Image.fromarray(base.astype(np.uint8), 'RGB')
img = img.filter(ImageFilter.GaussianBlur(radius=0.45))
img = img.resize((SIZE, SIZE), Image.Resampling.LANCZOS)
arr = np.asarray(img, dtype=np.float32)
fine = rng.normal(0.0, 2.2, size=arr.shape).astype(np.float32)
broad = gaussian_filter(rng.normal(size=(SIZE, SIZE)).astype(np.float32), sigma=22)
broad = broad[..., None]
arr += fine + 3.5 * broad / (np.std(broad) + 1e-6)
# Small interference-like colour tinting.
arr[..., 0] *= 1.012
arr[..., 1] *= 0.998
arr[..., 2] *= 1.014
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_v2.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()
