#!/usr/bin/env python3
"""Create a 4096x4096 allRGB image resembling fluorescent dopaminergic neurons.
Every RGB triplet occurs exactly once. Requires numpy, pillow and scipy.
Example: python dopamine_allrgb.py --seed 17 --output dopamine_allrgb.png
The scene is illustrative, not a real micrograph. Sorting the 24-bit palette
by luminance and matching it to the scene's luminance preserves structure while
ensuring the entire RGB gamut is used exactly once.
"""
import argparse
import gc
import numpy as np
from PIL import Image, ImageDraw, ImageFilter
from scipy.ndimage import gaussian_filter
SIZE = 4096
SMALL = 1024
N = SIZE * SIZE
def make_scene(seed):
rng = np.random.default_rng(seed)
# Draw at quarter resolution and supersample the finished luminous scene.
glow = Image.new('L', (SMALL, SMALL), 0)
sharp = Image.new('L', (SMALL, SMALL), 0)
gd, sd = ImageDraw.Draw(glow), ImageDraw.Draw(sharp)
# Curved axon/dendrite paths, with bead-like varicosities and soma shapes.
for _ in range(33):
x, y = rng.uniform(60, 970, 2)
theta = rng.uniform(0, 2 * np.pi)
length = int(rng.integers(80, 260))
points = [(x, y)]
for step in range(length):
theta += rng.normal(0, 0.085)
x += np.cos(theta) * rng.uniform(1.3, 2.8)
y += np.sin(theta) * rng.uniform(1.3, 2.8)
points.append((x, y))
if step % 12 == 0 and 0 <= x < SMALL and 0 <= y < SMALL:
r = rng.uniform(1.8, 5.0)
gd.ellipse((x-r*2, y-r*2, x+r*2, y+r*2), fill=205)
sd.ellipse((x-r, y-r, x+r, y+r), fill=240)
gd.line(points, fill=120, width=9, joint='curve')
sd.line(points, fill=195, width=2, joint='curve')
for _ in range(9):
x, y = rng.uniform(80, 940, 2)
r = rng.uniform(12, 31)
gd.ellipse((x-r*1.8, y-r*1.8, x+r*1.8, y+r*1.8), fill=220)
sd.ellipse((x-r, y-r, x+r, y+r), fill=225)
sd.ellipse((x-r*.42, y-r*.42, x+r*.42, y+r*.42), fill=45)
for _ in range(3600):
x, y = rng.uniform(0, SMALL, 2)
r = rng.uniform(.4, 2.5)
value = int(rng.integers(75, 230))
sd.ellipse((x-r, y-r, x+r, y+r), fill=value)
# Composite broad halo and crisp signal.
g = np.asarray(glow.filter(ImageFilter.GaussianBlur(13)), dtype=np.float32)
s = np.asarray(sharp.filter(ImageFilter.GaussianBlur(1.2)), dtype=np.float32)
field = np.clip(0.9*g + 1.05*s, 0, 255)
scene = Image.fromarray(field.astype(np.uint8), 'L').resize(
(SIZE, SIZE), Image.Resampling.BICUBIC)
del g, s, field
return np.asarray(scene, dtype=np.uint8).copy()
def render(seed, output):
rng = np.random.default_rng(seed + 1)
scene = make_scene(seed)
# A low-level uneven background plus film-like digital grain.
arr = scene.astype(np.float32)
del scene
# Dimming some bright processes preserves highlights in the allRGB remap.
arr *= 0.82
arr += rng.normal(0, 13, size=arr.shape).astype(np.float32)
arr += 16
np.clip(arr, 0, 255, out=arr)
target = arr.astype(np.uint8).ravel()
del arr
print('Sorting pixels by fluorescence intensity...', flush=True)
# Stable sort is unnecessary: randomized grain breaks up equal-valued regions.
pixel_order = np.argsort(target, kind='stable')
del target
gc.collect()
print('Sorting the complete RGB palette by perceived brightness...', flush=True)
palette_ids = np.arange(N, dtype=np.uint32)
r = (palette_ids >> 16).astype(np.uint16)
g = ((palette_ids >> 8) & 255).astype(np.uint16)
b = (palette_ids & 255).astype(np.uint16)
# Integer BT.709 brightness approximation. All RGB colors are retained.
luminosity = (54*r + 183*g + 19*b).astype(np.uint16)
del r, g, b
palette_order = np.argsort(luminosity, kind='stable')
del luminosity
gc.collect()
print('Distributing each RGB value exactly once...', flush=True)
# Keep peak RAM relatively low: uint32 packed pixels, then a uint8 RGB view.
packed = np.empty(N, dtype=np.uint32)
# Chunk to avoid an additional enormous temporary array.
chunk = 500_000
for start in range(0, N, chunk):
end = min(start + chunk, N)
packed[pixel_order[start:end]] = palette_ids[palette_order[start:end]]
del palette_ids, palette_order, pixel_order
gc.collect()
# Convert 0xRRGGBB packed numbers to an RGB byte array.
rgb = np.empty((N, 3), dtype=np.uint8)
rgb[:, 0] = packed >> 16
rgb[:, 1] = packed >> 8
rgb[:, 2] = packed
del packed
Image.fromarray(rgb.reshape(SIZE, SIZE, 3), 'RGB').save(output, optimize=False)
print(f'Saved {output}: {SIZE} x {SIZE}, all {N:,} RGB colors exactly once.')
def main():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--output', default='dopamine_allrgb.png')
parser.add_argument('--seed', type=int, default=17)
args = parser.parse_args()
render(args.seed, args.output)
if __name__ == '__main__':
main()
